WGCNA
Kristen Brown, modified ZD
1/24/2024
WGCNA Analysis
This is a Rmd file analyzing our raw count data by the WGCNA package as described in the manual and tutorials.
sessionInfo() #provides list of loaded packages and version of R. ## R version 4.3.2 (2023-10-31)
## Platform: aarch64-apple-darwin20 (64-bit)
## Running under: macOS Ventura 13.0
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.11.0
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: America/New_York
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## loaded via a namespace (and not attached):
## [1] digest_0.6.35 R6_2.5.1 fastmap_1.2.0 xfun_0.43
## [5] cachem_1.0.8 knitr_1.45 htmltools_0.5.8.1 rmarkdown_2.26
## [9] lifecycle_1.0.4 cli_3.6.2 sass_0.4.9 jquerylib_0.1.4
## [13] compiler_4.3.2 rstudioapi_0.16.0 tools_4.3.2 evaluate_0.23
## [17] bslib_0.7.0 yaml_2.3.8 rlang_1.1.3 jsonlite_1.8.8First, download the necessary packages.
install.packages("BiocManager")
library("BiocManager")
BiocManager::install("impute", type = "source")
BiocManager::install("WGCNA",force = TRUE)
BiocManager::install("vsn")
install.packages("dendextend")Next, load the packages
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(tidyr)
library(Rmisc)## Loading required package: lattice## Loading required package: plyr## ------------------------------------------------------------------------------## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)## ------------------------------------------------------------------------------##
## Attaching package: 'plyr'## The following objects are masked from 'package:dplyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarizelibrary(lubridate)##
## Attaching package: 'lubridate'## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, unionlibrary(ggplot2)
library(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recodelibrary(MASS)##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## selectlibrary(mgcv)## Loading required package: nlme##
## Attaching package: 'nlme'## The following object is masked from 'package:dplyr':
##
## collapse## This is mgcv 1.9-1. For overview type 'help("mgcv-package")'.library(MuMIn)
library(emmeans)
library(locfit)## locfit 1.5-9.9 2024-03-01library(DESeq2)## Loading required package: S4Vectors## Loading required package: stats4## Loading required package: BiocGenerics##
## Attaching package: 'BiocGenerics'## The following objects are masked from 'package:lubridate':
##
## intersect, setdiff, union## The following objects are masked from 'package:dplyr':
##
## combine, intersect, setdiff, union## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs## The following objects are masked from 'package:base':
##
## anyDuplicated, aperm, append, as.data.frame, basename, cbind,
## colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
## get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
## match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
## Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort,
## table, tapply, union, unique, unsplit, which.max, which.min##
## Attaching package: 'S4Vectors'## The following objects are masked from 'package:lubridate':
##
## second, second<-## The following object is masked from 'package:plyr':
##
## rename## The following object is masked from 'package:tidyr':
##
## expand## The following objects are masked from 'package:dplyr':
##
## first, rename## The following object is masked from 'package:utils':
##
## findMatches## The following objects are masked from 'package:base':
##
## expand.grid, I, unname## Loading required package: IRanges##
## Attaching package: 'IRanges'## The following object is masked from 'package:nlme':
##
## collapse## The following object is masked from 'package:lubridate':
##
## %within%## The following object is masked from 'package:plyr':
##
## desc## The following objects are masked from 'package:dplyr':
##
## collapse, desc, slice## Loading required package: GenomicRanges## Loading required package: GenomeInfoDb## Loading required package: SummarizedExperiment## Loading required package: MatrixGenerics## Loading required package: matrixStats##
## Attaching package: 'matrixStats'## The following object is masked from 'package:plyr':
##
## count## The following object is masked from 'package:dplyr':
##
## count##
## Attaching package: 'MatrixGenerics'## The following objects are masked from 'package:matrixStats':
##
## colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
## colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
## colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
## colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
## colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
## colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
## colWeightedMeans, colWeightedMedians, colWeightedSds,
## colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
## rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
## rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
## rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
## rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
## rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
## rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
## rowWeightedSds, rowWeightedVars## Loading required package: Biobase## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.##
## Attaching package: 'Biobase'## The following object is masked from 'package:MatrixGenerics':
##
## rowMedians## The following objects are masked from 'package:matrixStats':
##
## anyMissing, rowMedianslibrary(impute)
library(WGCNA)## Loading required package: dynamicTreeCut## Loading required package: fastcluster##
## Attaching package: 'fastcluster'## The following object is masked from 'package:stats':
##
## hclust## ##
## Attaching package: 'WGCNA'## The following object is masked from 'package:IRanges':
##
## cor## The following object is masked from 'package:S4Vectors':
##
## cor## The following object is masked from 'package:stats':
##
## corlibrary(xfun)##
## Attaching package: 'xfun'## The following objects are masked from 'package:base':
##
## attr, isFALSElibrary(genefilter)##
## Attaching package: 'genefilter'## The following objects are masked from 'package:MatrixGenerics':
##
## rowSds, rowVars## The following objects are masked from 'package:matrixStats':
##
## rowSds, rowVars## The following object is masked from 'package:MASS':
##
## area## The following object is masked from 'package:car':
##
## Anovalibrary(vsn)
library(RColorBrewer)
library(pheatmap)
library(dendextend)##
## ---------------------
## Welcome to dendextend version 1.17.1
## Type citation('dendextend') for how to cite the package.
##
## Type browseVignettes(package = 'dendextend') for the package vignette.
## The github page is: https://github.com/talgalili/dendextend/
##
## Suggestions and bug-reports can be submitted at: https://github.com/talgalili/dendextend/issues
## You may ask questions at stackoverflow, use the r and dendextend tags:
## https://stackoverflow.com/questions/tagged/dendextend
##
## To suppress this message use: suppressPackageStartupMessages(library(dendextend))
## ---------------------##
## Attaching package: 'dendextend'## The following object is masked from 'package:stats':
##
## cutreesessionInfo()## R version 4.3.2 (2023-10-31)
## Platform: aarch64-apple-darwin20 (64-bit)
## Running under: macOS Ventura 13.0
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.11.0
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: America/New_York
## tzcode source: internal
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] dendextend_1.17.1 pheatmap_1.0.12
## [3] RColorBrewer_1.1-3 vsn_3.68.0
## [5] genefilter_1.82.1 xfun_0.43
## [7] WGCNA_1.72-5 fastcluster_1.2.6
## [9] dynamicTreeCut_1.63-1 impute_1.74.1
## [11] DESeq2_1.40.2 SummarizedExperiment_1.30.2
## [13] Biobase_2.60.0 MatrixGenerics_1.12.3
## [15] matrixStats_1.3.0 GenomicRanges_1.52.1
## [17] GenomeInfoDb_1.36.4 IRanges_2.34.1
## [19] S4Vectors_0.38.2 BiocGenerics_0.46.0
## [21] locfit_1.5-9.9 emmeans_1.10.1
## [23] MuMIn_1.47.5 mgcv_1.9-1
## [25] nlme_3.1-164 MASS_7.3-60.0.1
## [27] car_3.1-2 carData_3.0-5
## [29] ggplot2_3.5.1 lubridate_1.9.3
## [31] Rmisc_1.5.1 plyr_1.8.9
## [33] lattice_0.22-6 tidyr_1.3.1
## [35] dplyr_1.1.4
##
## loaded via a namespace (and not attached):
## [1] rstudioapi_0.16.0 jsonlite_1.8.8 magrittr_2.0.3
## [4] TH.data_1.1-2 estimability_1.5.1 rmarkdown_2.26
## [7] zlibbioc_1.46.0 vctrs_0.6.5 memoise_2.0.1
## [10] RCurl_1.98-1.14 base64enc_0.1-3 htmltools_0.5.8.1
## [13] S4Arrays_1.0.6 Formula_1.2-5 sass_0.4.9
## [16] bslib_0.7.0 htmlwidgets_1.6.4 sandwich_3.1-0
## [19] zoo_1.8-12 cachem_1.0.8 lifecycle_1.0.4
## [22] iterators_1.0.14 pkgconfig_2.0.3 Matrix_1.6-5
## [25] R6_2.5.1 fastmap_1.2.0 GenomeInfoDbData_1.2.10
## [28] digest_0.6.35 colorspace_2.1-0 AnnotationDbi_1.62.2
## [31] Hmisc_5.1-2 RSQLite_2.3.6 fansi_1.0.6
## [34] timechange_0.3.0 httr_1.4.7 abind_1.4-5
## [37] compiler_4.3.2 bit64_4.0.5 withr_3.0.0
## [40] doParallel_1.0.17 htmlTable_2.4.2 backports_1.5.0
## [43] BiocParallel_1.34.2 viridis_0.6.5 DBI_1.2.2
## [46] DelayedArray_0.26.7 tools_4.3.2 foreign_0.8-86
## [49] nnet_7.3-19 glue_1.7.0 grid_4.3.2
## [52] checkmate_2.3.1 cluster_2.1.6 generics_0.1.3
## [55] gtable_0.3.5 preprocessCore_1.62.1 data.table_1.15.4
## [58] utf8_1.2.4 XVector_0.40.0 foreach_1.5.2
## [61] pillar_1.9.0 stringr_1.5.1 limma_3.56.2
## [64] splines_4.3.2 survival_3.6-4 bit_4.0.5
## [67] annotate_1.78.0 tidyselect_1.2.1 GO.db_3.17.0
## [70] Biostrings_2.68.1 knitr_1.45 gridExtra_2.3
## [73] stringi_1.8.3 yaml_2.3.8 evaluate_0.23
## [76] codetools_0.2-20 tibble_3.2.1 BiocManager_1.30.23
## [79] cli_3.6.2 affyio_1.70.0 rpart_4.1.23
## [82] xtable_1.8-4 munsell_0.5.1 jquerylib_0.1.4
## [85] Rcpp_1.0.12 coda_0.19-4.1 png_0.1-8
## [88] XML_3.99-0.16.1 parallel_4.3.2 blob_1.2.4
## [91] bitops_1.0-7 viridisLite_0.4.2 mvtnorm_1.2-5
## [94] scales_1.3.0 affy_1.78.2 purrr_1.0.2
## [97] crayon_1.5.2 rlang_1.1.3 KEGGREST_1.40.1
## [100] multcomp_1.4-25Load in raw count data
cts_raw <- read.csv("../../TagSeq_Output/HeronPdam_gene_count_matrix.csv") #load in data
rownames(cts_raw) <- cts_raw[,1] #set first column that contains gene names as rownames
cts_raw <- cts_raw[,-1] #remove the column with gene names
head(cts_raw)## RF13B_S85_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1 3
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1 2
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1 0
## Pocillopora_acuta_HIv2___TS.g3788.t1 0
## RF13D_S112_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1 8
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1 0
## Pocillopora_acuta_HIv2___TS.g3788.t1 0
## RF14B_S116_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1 4
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1 0
## Pocillopora_acuta_HIv2___TS.g3788.t1 0
## RF14C_S95_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1 0
## Pocillopora_acuta_HIv2___TS.g3788.t1 1
## RF15B_S122_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1 2
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1 0
## Pocillopora_acuta_HIv2___TS.g3788.t1 0
## RF15D_S84_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1 0
## Pocillopora_acuta_HIv2___TS.g3788.t1 0
## RF16A_S118_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1 0
## Pocillopora_acuta_HIv2___TS.g3788.t1 0
## RF16C_S97_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1 0
## Pocillopora_acuta_HIv2___TS.g3788.t1 0
## RF17B_S80_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1 10
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1 0
## Pocillopora_acuta_HIv2___TS.g3788.t1 1
## RF17D_S99_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1 0
## Pocillopora_acuta_HIv2___TS.g3788.t1 1
## RF18B_S90_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1 0
## Pocillopora_acuta_HIv2___TS.g3788.t1 0
## RF18D_S113_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1 0
## Pocillopora_acuta_HIv2___TS.g3788.t1 0
## RF19B_S103_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1 0
## Pocillopora_acuta_HIv2___TS.g3788.t1 1
## RF19C_S83_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1 0
## Pocillopora_acuta_HIv2___TS.g3788.t1 0
## RF20B_S88_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1 2
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1 0
## Pocillopora_acuta_HIv2___TS.g3788.t1 0
## RF20C_S115_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1 0
## Pocillopora_acuta_HIv2___TS.g3788.t1 1
## RF22B_S91_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1 0
## Pocillopora_acuta_HIv2___TS.g3788.t1 1
## RF22C_S126_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1 0
## Pocillopora_acuta_HIv2___TS.g3788.t1 1
## RF23A_S125_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1 5
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1 6
## Pocillopora_acuta_HIv2___TS.g3788.t1 0
## RF23C_S108_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1 0
## Pocillopora_acuta_HIv2___TS.g3788.t1 0
## RF24B_S92_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1 0
## Pocillopora_acuta_HIv2___TS.g3788.t1 0
## RF24D_S111_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1 7
## Pocillopora_acuta_HIv2___TS.g3788.t1 0
## RF25A_S96_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1 0
## Pocillopora_acuta_HIv2___TS.g3788.t1 1
## RF25C_S110_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1 0
## Pocillopora_acuta_HIv2___TS.g3788.t1 1
## RS11B_S86_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1 1
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1 6
## Pocillopora_acuta_HIv2___TS.g3788.t1 1
## RS11D_S87_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1 0
## Pocillopora_acuta_HIv2___TS.g3788.t1 0
## RS12A_S105_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1 0
## Pocillopora_acuta_HIv2___TS.g3788.t1 1
## RS12C_S106_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1 2
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1 0
## Pocillopora_acuta_HIv2___TS.g3788.t1 0
## RS13A_S121_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1 0
## Pocillopora_acuta_HIv2___TS.g3788.t1 0
## RS13C_S119_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1 2
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1 0
## Pocillopora_acuta_HIv2___TS.g3788.t1 1
## RS14B_S107_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1 0
## Pocillopora_acuta_HIv2___TS.g3788.t1 0
## RS14C_S109_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1 8
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1 0
## Pocillopora_acuta_HIv2___TS.g3788.t1 0
## RS15B_S104_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1 8
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1 0
## Pocillopora_acuta_HIv2___TS.g3788.t1 0
## RS15D_S82_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1 10
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1 0
## Pocillopora_acuta_HIv2___TS.g3788.t1 0
## RS1B_S102_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1 1
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1 1
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1 0
## Pocillopora_acuta_HIv2___TS.g3788.t1 0
## RS1C_S127_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1 1
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1 4
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1 0
## Pocillopora_acuta_HIv2___TS.g3788.t1 1
## RS2B_S123_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1 0
## Pocillopora_acuta_HIv2___TS.g3788.t1 0
## RS2C_S89_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1 0
## Pocillopora_acuta_HIv2___TS.g3788.t1 1
## RS3B_S100_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1 0
## Pocillopora_acuta_HIv2___TS.g3788.t1 0
## RS3D_S117_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1 0
## Pocillopora_acuta_HIv2___TS.g3788.t1 0
## RS6A_S81_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1 0
## Pocillopora_acuta_HIv2___TS.g3788.t1 0
## RS6D_S124_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1 2
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1 0
## Pocillopora_acuta_HIv2___TS.g3788.t1 0
## RS7B_S120_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1 0
## Pocillopora_acuta_HIv2___TS.g3788.t1 0
## RS7C_S101_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1 0
## Pocillopora_acuta_HIv2___TS.g3788.t1 0
## RS8B_S94_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1 5
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1 0
## Pocillopora_acuta_HIv2___TS.g3788.t1 0
## RS8C_S114_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1 0
## Pocillopora_acuta_HIv2___TS.g3788.t1 0
## RS9A_S98_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1 10
## Pocillopora_acuta_HIv2___TS.g3788.t1 1
## RS9C_S93_ALL.bam.gtf
## Pocillopora_acuta_HIv2___RNAseq.g23616.t1 2
## Pocillopora_acuta_HIv2___RNAseq.g4197.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g12137.t1 8
## Pocillopora_acuta_HIv2___RNAseq.g4145.t1 0
## Pocillopora_acuta_HIv2___RNAseq.g4784.t1 5
## Pocillopora_acuta_HIv2___TS.g3788.t1 1Clean up sample names from “RF13B_S85_ALL.bam.gtf” to “RF13B”
colnames(cts_raw) <- gsub("_[^_]+$", "",colnames(cts_raw)) #get rid of "_ALL.bam.gtf"
colnames(cts_raw) <- gsub("_[^_]+$", "",colnames(cts_raw)) #get rid of "_S85"
head(colnames(cts_raw)) #see first 6 clean sample names## [1] "RF13B" "RF13D" "RF14B" "RF14C" "RF15B" "RF15D"Metadata from this dataset
coldata <- read.csv("../../../TagSeq_Submission/RNA Submission Sample List metadata.csv") #read in metadata file
coldata <- plyr::rename(coldata, c("Sample.Name"="Coral_ID")) #Make a column that represents the colonies
coldata$Colony <- gsub("A", "", coldata$Coral_ID) #ID which sample is from which colony by removing letters after colony name
coldata$Colony <- gsub("B", "", coldata$Colony) #ID which sample is from which colony by removing letters after colony name
coldata$Colony <- gsub("C", "", coldata$Colony) #ID which sample is from which colony by removing letters after colony name
coldata$Colony <- gsub("D", "", coldata$Colony) #ID which sample is from which colony by removing letters after colony name
coldata$Origin <- factor(coldata$Origin, levels = c("Slope","Flat")) #set variables to factors, with Slope as the baseline
coldata$Colony <- factor(coldata$Colony) #set variables to factors
coldata$Treatment <- factor(coldata$Treatment) #set variables to factors
head(coldata)## Coral_ID Origin Treatment Colony
## 1 RF13B Flat Variable RF13
## 2 RF13D Flat Stable RF13
## 3 RF14B Flat Variable RF14
## 4 RF14C Flat Stable RF14
## 5 RF15B Flat Variable RF15
## 6 RF15D Flat Stable RF15Data sanity checks:
all(rownames(coldata$Coral_ID) %in% colnames(cts_raw)) #are all of the sample names (rows) in the metadata column names in the gene count matrix?## [1] TRUEall(rownames(coldata$Coral_ID) == colnames(cts_raw)) #are they the same in the same order?## [1] TRUEFiltering
dim(cts_raw)## [1] 33730 48cts_filt <- cts_raw[rowSums(!as.matrix(cts_raw)) < ncol(cts_raw), ] #here we remove all genes that were not expressed in any of our samples- we are left with 24,233 genes.
dim(cts_filt)## [1] 24233 48pOverA filtering to reduce dataset
ffun <- filterfun(pOverA(0.25,10)) #set up filtering parameters
cts_filtpoa <- genefilter((cts_filt), ffun) #apply filter
sum(cts_filtpoa) #count number of genes left## [1] 9056cts_filtpoa <- cts_filt[cts_filtpoa,] #keep only rows that passed filter, 9056 genes
dim(cts_filtpoa)## [1] 9056 48library("DESeq2")
dds <- DESeqDataSetFromMatrix(countData = cts_filtpoa,
colData=coldata,
design= ~Origin+Treatment)
dds## class: DESeqDataSet
## dim: 9056 48
## metadata(1): version
## assays(1): counts
## rownames(9056): Pocillopora_acuta_HIv2___RNAseq.g27841.t1
## Pocillopora_acuta_HIv2___RNAseq.g14011.t1 ...
## Pocillopora_acuta_HIv2___RNAseq.g27413.t1
## Pocillopora_acuta_HIv2___RNAseq.g10545.t1
## rowData names(0):
## colnames(48): RF13B RF13D ... RS9A RS9C
## colData names(4): Coral_ID Origin Treatment Colony# transform the data
vst <- vst(dds, blind=FALSE) # transform it vstEffects of transformations on the variance
meanSdPlot(assay(vst))
Data quality assessment by sample clustering and visualization
Heatmap of the sample-to-sample distances
sampleDists <- dist(t(assay(vst)))sampleDistMatrix <- as.matrix(sampleDists)
rownames(sampleDistMatrix) <- paste(vst$Origin, vst$Treatment, vst$Coral_ID, sep="-")
colnames(sampleDistMatrix) <- NULL
colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)
pheatmap(sampleDistMatrix,
clustering_distance_rows=sampleDists,
clustering_distance_cols=sampleDists,
col=colors)
Principal component plot of the samples
withoutliers <- plotPCA(vst, intgroup=c("Origin", "Treatment"))
withoutliers
Removing outlier “RF16A” and “RF16C”
cts_raw_outrm <- cts_raw %>% dplyr::select(-c("RF16A","RF16C")) #remove those columns from count matrix
coldata_outrm <- coldata %>% filter(Coral_ID != c("RF16A","RF16C")) #removed those rows from metadataUsing filtered data instead of raw counts as input
cts_filt_outrm <- cts_raw_outrm[rowSums(!as.matrix(cts_raw_outrm)) < ncol(cts_raw_outrm), ] #here we remove all genes that were not expressed in any of our samples- we are left with 24,220 genes.
dim(cts_filt_outrm)## [1] 24220 46pOverA filtering to reduce dataset
ffun <- filterfun(pOverA(0.25,10)) #set up filtering parameters
cts_filt_outrm_poa <- genefilter((cts_filt_outrm), ffun) #apply filter
sum(cts_filt_outrm_poa) #count number of genes left## [1] 9012cts_filt_outrm_poa <- cts_filt_outrm[cts_filt_outrm_poa,] #keep only rows that passed filter, 9011 genes
cts_filtpoa <- cts_filt_outrm_poaData sanity checks:
all(rownames(coldata_outrm$Coral_ID) %in% colnames(cts_filt_outrm_poa)) #are all of the sample names (rows) in the metadata column names in the gene count matrix?## [1] TRUEall(rownames(coldata_outrm$Coral_ID) == colnames(cts_filt_outrm_poa)) #are they the same in the same order?## [1] TRUEwrite.csv(cts_filt_outrm_poa, "../../output/Filtered_gene_count_matrix.csv")library("DESeq2")
dds2 <- DESeqDataSetFromMatrix(countData = cts_filtpoa,
colData=coldata_outrm,
design= ~Origin+Treatment)
dds2## class: DESeqDataSet
## dim: 9012 46
## metadata(1): version
## assays(1): counts
## rownames(9012): Pocillopora_acuta_HIv2___RNAseq.g27841.t1
## Pocillopora_acuta_HIv2___RNAseq.g14011.t1 ...
## Pocillopora_acuta_HIv2___RNAseq.g27413.t1
## Pocillopora_acuta_HIv2___RNAseq.g10545.t1
## rowData names(0):
## colnames(46): RF13B RF13D ... RS9A RS9C
## colData names(4): Coral_ID Origin Treatment Colony# transform the data
vst2 <- vst(dds2, blind=FALSE) # transform it vstEffects of transformations on the variance
meanSdPlot(assay(vst2))
Data quality assessment by sample clustering and visualization
Heatmap of the sample-to-sample distances
sampleDists <- dist(t(assay(vst2)))library("RColorBrewer")
sampleDistMatrix <- as.matrix(sampleDists)
rownames(sampleDistMatrix) <- paste(vst2$Origin, vst2$Treatment, vst2$Coral_ID, sep="-")
colnames(sampleDistMatrix) <- NULL
colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)
pheatmap(sampleDistMatrix,
clustering_distance_rows=sampleDists,
clustering_distance_cols=sampleDists,
col=colors)
###Principal component plot of the samples
withoutoutliers <- plotPCA(vst2, intgroup=c("Origin", "Treatment"))
withoutoutliers
compare_figs <- cowplot::plot_grid(withoutliers, withoutoutliers,ncol=2, align="h")
compare_figs
Transpose the filtered gene count matrix so that the gene IDs are rows and the sample IDs are columns.
# transform the data
vst2 <- vst(dds2) # transform it vst
vst2 <- assay(vst2) # call only the transformed coutns in the dds object
vst2 <- t(vst2) # transpose columns to rows and vice versaCheck for genes and samples with too many missing values with goodSamplesGenes. There shouldn’t be any because we performed pre-filtering
dim(vst2) # 9012 genes; 46 samples## [1] 46 9012gsg <- goodSamplesGenes(vst2, verbose = 3); # We first check for genes and samples with too many missing values## Flagging genes and samples with too many missing values...
## ..step 1gsg$allOK # If the statement returns TRUE, all genes have passed the cuts. ## [1] TRUEConduct WGCNA
library(WGCNA)
# The following setting is important, do not omit.
options(stringsAsFactors = FALSE)##Soft threshold
dim(vst2) # 46 9012## [1] 46 9012# Choose a set of soft-thresholding powers
powers = c(c(1:10), seq(from = 12, to=20, by=2))
# Call the network topology analysis function
#the below takes a long time to run, so is commented out and the pre-run results are loaded in below
## sft = pickSoftThreshold(vst2, powerVector = powers, verbose = 5) #...wait for this to finish
## save(sft, file = "../../output/WGCNA/sft.RData")
load("../../output/WGCNA/sft.RData")
# pickSoftThreshold
# performs the analysis of network topology and aids the
# user in choosing a proper soft-thresholding power.
# The user chooses a set of candidate powers (the function provides suitable default values)
# function returns a set of network indices that should be inspected
sizeGrWindow(9, 5) # set window size
# png to output
png("../../output/WGCNA/sft.png", 1000, 1000, pointsize=20)
par(mfrow = c(1,2));
cex1 = 0.9;
# Scale-free topology fit index as a function of the soft-thresholding power
plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n",
main = paste("Scale independence"));
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
labels=powers,cex=cex1,col="red");
# this line corresponds to using an R^2 cut-off of h
abline(h=0.90,col="red")
# Mean connectivity as a function of the soft-thresholding power
plot(sft$fitIndices[,1], sft$fitIndices[,5],
xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
main = paste("Mean connectivity"))
text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red")
dev.off() # output ## pdf
## 2#I used a scale-free topology fit index **R^2 of 0.9**. This lowest recommended R^2 by Langfelder and Horvath is 0.8. I chose 0.9 because we want to use the smallest soft thresholding power that maximizes with model fit. It appears that our **soft thresholding power is 5** because it is the lowest power above the R^2=0.9 threshold that maximizes with model fit. Look for outliers by examining tree of samples
sampleTree = hclust(dist(vst2), method = "average") # Next we cluster the samples (in contrast to clustering genes that will come later) to see if there are any obvious outliers.There don't look to be any outliers, so we will move on with business as usual.
sizeGrWindow(12,9)
par(cex = 0.6);
par(mar = c(0,4,2,0))
plot(sampleTree)Start the step-wise module construction:
Step 1: Create adjacency matrix
softPower = 5 # set the soft threshold based on the plots above
# signed
#to get this to work with a ton of genes >30K, I had to increase my memory limit using: memory.limit(size = 35000)
adjacency_sign = adjacency(vst2, power = softPower, type="signed") #Calculate adjacencyStep 2: Turn adjacency into topological overlap: Calculation of the topological overlap matrix, (TOM) and the corresponding dissimilarity, from a given adjacency matrix.
#the below takes a long time to run, so is commented out and the pre-run results are loaded in below
## TOM_sign = TOMsimilarity(adjacency_sign, TOMType="signed") #Translate adjacency into topological overlap matrix
## save(TOM_sign, file = "../../output/WGCNA/TOM_sign.Rdata")
load("../../output/WGCNA/TOM_sign.Rdata")
dissTOM_sign = 1-TOM_signStep 3: Call the hierarchical clustering function - plot the tree
# Call the hierarchical clustering function
#to get this to work, I had to increase my memory limit using: memory.limit(size = 45000)
geneTree_sign = hclust(as.dist(dissTOM_sign), method = "average");
# Plot the resulting clustering tree (dendrogram) Each leaf corresponds to a gene, branches grouping together densely are interconnected, highly co-expressed genes.
sizeGrWindow(12,9)
plot(geneTree_sign, xlab="", sub="", main = "Gene clustering on TOM-based dissimilarity - SIGNED",
labels = FALSE, hang = 0.04);Step 4: Set module size and ‘cutreeDynamic’ to create clusters
#Module identification is essentially cutting the branches off the tree in the dendrogram above. We like large modules, so we set the **minimum module size** relatively high, so we will set the minimum size at 30.
minModuleSize = 30; # set this for the subseqent call...
dynamicMods_sign = cutreeDynamic(dendro = geneTree_sign, distM = dissTOM_sign,
deepSplit = 1, pamRespectsDendro = FALSE,
minClusterSize = minModuleSize);## ..cutHeight not given, setting it to 0.959 ===> 99% of the (truncated) height range in dendro.
## ..done.table(dynamicMods_sign) # number of genes per module. Module 0 is reserved for unassigned genes. The are other modules will be listed largest to smallest. ## dynamicMods_sign
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
## 2558 1809 942 917 500 425 396 220 219 218 180 170 154 105 71 65
## 17
## 63Step 5: convert numeric network to colors and plot the dendrogram
# Convert numeric lables into colors
dynamicColors_sign = labels2colors(dynamicMods_sign) # add colors to module labels (previously numbers)
table(dynamicColors_sign) # lets look at this table...## dynamicColors_sign
## black blue brown cyan green greenyellow
## 396 1809 942 105 500 180
## grey60 lightcyan magenta midnightblue pink purple
## 63 65 219 71 220 218
## red salmon tan turquoise yellow
## 425 154 170 2558 917# Plot the dendrogram and colors underneath
plotDendroAndColors(geneTree_sign, dynamicColors_sign, "Dynamic Tree Cut - SIGNED",
dendroLabels = FALSE, hang = 0.03,
addGuide = TRUE, guideHang = 0.05,
main = "Gene dendrogram and module colors 'SIGNED'")
sizeGrWindow(8,6)
png("../../output/WGCNA/GeneDendrogram.png", 1000, 1000, pointsize=20)
plotDendroAndColors(geneTree_sign, dynamicColors_sign, "Dynamic Tree Cut - SIGNED",
dendroLabels = FALSE, hang = 0.03,
addGuide = TRUE, guideHang = 0.05,
main = "Gene dendrogram and module colors 'SIGNED'")
dev.off()## pdf
## 2Step 6: Calculate Eigengenes - view thier connectivity based on ‘MEDiss = 1-cor(MEs)’
# Calculate eigengenes
MEList = moduleEigengenes(vst2, colors = dynamicColors_sign, softPower = 5)
MEs = MEList$eigengenes
MEs## MEblack MEblue MEbrown MEcyan MEgreen
## RF13B 0.035738890 -0.0005614067 0.086214000 0.0230457125 0.06457070
## RF13D 0.016715232 0.0355808790 0.077795819 0.0108599933 0.04977675
## RF14B -0.146566530 -0.0989830172 0.243945678 0.0089738153 0.09690707
## RF14C -0.144417844 0.0293510325 0.231391380 0.0555552483 0.12492062
## RF15B 0.361363331 0.0125541542 0.195194429 0.0669848895 0.09848360
## RF15D 0.340691436 -0.0628815348 0.130600265 0.0613695551 0.15041068
## RF17B -0.166389167 -0.0083085168 0.198989474 -0.0330477836 0.06580995
## RF17D -0.174336507 0.0274499964 0.242761454 0.0109581641 0.20156729
## RF18B -0.163062969 0.0730460426 0.213165173 0.0428568437 0.15059016
## RF18D -0.198796035 0.1501095732 0.221598881 0.0156327910 0.26172808
## RF19B -0.052085383 0.1809702199 0.156245872 0.0534562058 0.11871425
## RF19C -0.133192276 -0.1908284494 0.207082710 0.0058409291 0.09057738
## RF20B -0.203191630 0.1418285398 0.167736368 0.0175489549 0.09355805
## RF20C -0.079445376 0.0589582967 0.220340540 0.0517627628 0.08359879
## RF22B 0.013623481 -0.1171079552 0.003985773 0.0146616850 -0.06040064
## RF22C -0.011991191 0.1497873255 0.042606148 0.0054811614 -0.06439484
## RF23A 0.312454266 0.1092940451 0.061065835 0.0347675690 0.09520240
## RF23C 0.295349869 0.0345798457 0.114909745 0.0202737297 0.03245128
## RF24B 0.356715054 0.1929708845 0.106044589 0.0913070182 0.03457393
## RF24D 0.306302794 0.1516425029 0.076269408 0.0729544308 0.08671047
## RF25A 0.080575924 -0.0834460336 -0.007865851 -0.0080197431 -0.04258953
## RF25C 0.093458065 0.1069453692 -0.019058548 0.0240384357 0.09631201
## RS11B -0.004644330 -0.1099408605 -0.092245924 0.0376314369 0.10628308
## RS11D -0.026519519 0.0766325673 -0.153945434 0.0009392144 0.34571416
## RS12A -0.018592944 0.0934519345 -0.215804220 0.0079270236 -0.25126870
## RS12C 0.025448390 -0.2497580080 -0.180395235 0.0528064649 -0.06947225
## RS13A -0.013832874 0.2242246956 -0.087780763 0.0043265989 -0.17859214
## RS13C -0.016960538 0.0901775914 -0.092445893 0.0264009339 -0.20267882
## RS14B -0.025815615 -0.1101678607 -0.138277199 0.0152473895 -0.14673511
## RS14C -0.093490223 -0.3533046553 -0.224348383 -0.0853239910 -0.29804769
## RS15B -0.043303446 -0.0724784483 -0.135120081 0.0254364294 -0.25047579
## RS15D -0.063247364 -0.4464606863 -0.250028576 0.0434992663 -0.21085363
## RS1B -0.071734949 0.1666814228 -0.146538808 0.0318043922 0.01154907
## RS1C -0.081189863 -0.0208045578 -0.082597955 0.0382703168 0.03230402
## RS2B 0.081141470 0.0226962013 -0.046649316 -0.0191289688 -0.10542840
## RS2C 0.001196774 0.0341686249 -0.061856762 0.0202354959 -0.13064096
## RS3B -0.007344414 0.0096703526 -0.121890691 0.0118973832 -0.01743041
## RS3D -0.005781953 0.0046080335 -0.079855011 0.0442362556 0.03886139
## RS6A -0.070105162 -0.0231943131 -0.147206872 -0.9641101677 -0.14047140
## RS6D -0.055100416 -0.1924917784 -0.048878771 -0.0740728504 -0.01160300
## RS7B -0.115318529 -0.0362586532 -0.189674392 -0.0186040960 -0.26884057
## RS7C -0.191763564 -0.3670929579 -0.153268154 0.0157100485 -0.22995477
## RS8B 0.035698058 0.1335243682 -0.069128256 0.0594849913 -0.01327168
## RS8C 0.017104694 0.1220210827 -0.089764315 0.0442748956 -0.03140854
## RS9A -0.008609524 0.0942830753 -0.113018368 0.0003933584 0.23592747
## RS9C 0.013252409 0.0168610358 -0.050299764 0.0334558101 -0.04254375
## MEgreenyellow MEgrey60 MElightcyan MEmagenta MEmidnightblue
## RF13B 0.03096313 0.048813655 0.004539119 0.040533178 -0.0245491734
## RF13D 0.07466948 -0.002747228 -0.048921622 0.004749849 -0.0183763055
## RF14B -0.02351578 0.052054548 0.007711938 -0.055936540 -0.0228806751
## RF14C 0.05431928 0.103521196 -0.069472828 -0.015987095 0.0063273755
## RF15B 0.11598069 -0.027006621 -0.344930466 -0.362117043 0.1513293329
## RF15D -0.03629918 0.005716079 -0.167570877 -0.295093407 0.2593929751
## RF17B 0.04043036 0.037559058 -0.061350955 -0.105576226 0.0522009454
## RF17D 0.08347539 0.240079479 -0.003870975 -0.071316708 0.0889135020
## RF18B 0.17271940 0.088913590 -0.120851396 -0.157907924 0.0770626674
## RF18D 0.16454055 0.312762543 -0.045673177 -0.139895717 0.1822861472
## RF19B 0.11959708 0.096034394 0.029738991 -0.021384698 -0.0659116912
## RF19C -0.06142188 0.035342775 -0.226848538 -0.113708482 0.0063878362
## RF20B 0.20131380 0.085070425 -0.066774285 -0.071679981 -0.0725077225
## RF20C 0.03843045 0.065751826 -0.005878656 -0.079070138 -0.0267314858
## RF22B -0.12513476 -0.112608524 0.034288078 -0.033064166 -0.0109426095
## RF22C 0.18911651 -0.134875929 0.010881417 -0.037849207 -0.1930922669
## RF23A 0.07443417 -0.015840141 -0.080743552 -0.255329622 0.1073787126
## RF23C 0.12020491 -0.135276896 -0.256935676 -0.350307537 0.0506762137
## RF24B 0.16937663 -0.090056137 -0.119388795 -0.241229366 0.0014127112
## RF24D 0.11689102 0.002523370 0.002249586 -0.107578383 0.0916879301
## RF25A -0.11865150 -0.047263339 0.136314599 0.060319789 0.0447707556
## RF25C 0.06960809 0.065371487 0.115822838 0.091429118 0.1221668700
## RS11B -0.12420996 0.061302028 -0.150987001 0.022147377 0.1402703642
## RS11D 0.05539636 0.541732945 -0.060301638 -0.009668998 0.4557546749
## RS12A -0.07207071 -0.112117851 0.410469539 0.213579019 -0.2390647714
## RS12C -0.29590257 -0.015964786 0.110170934 0.097511706 0.0886475060
## RS13A 0.13529317 -0.157437914 0.243742958 0.153568478 -0.3479725582
## RS13C 0.05040341 -0.225760997 0.136072703 0.120176386 -0.2901638771
## RS14B -0.17606312 -0.134924778 0.149095999 0.126931780 -0.1067129932
## RS14C -0.32994363 -0.250860431 0.111490397 0.110626775 -0.1183388272
## RS15B -0.10615944 -0.212941292 0.118437262 0.196764982 -0.2220598840
## RS15D -0.50142342 -0.122195112 -0.002867751 0.192810781 0.0205032779
## RS1B 0.09307705 0.102779430 0.092266234 0.165333200 -0.0786940195
## RS1C 0.03100818 0.073828989 -0.122813217 0.021896321 0.1067984725
## RS2B 0.03704151 -0.153620437 -0.026877271 -0.043371606 -0.0752202101
## RS2C 0.02758929 -0.078462245 0.011979459 0.049463264 0.0101279467
## RS3B 0.02123936 -0.033816965 -0.028347890 -0.084532348 -0.0601723773
## RS3D 0.01938915 0.033029583 -0.081046733 -0.004556133 -0.0001098878
## RS6A 0.07246209 -0.073104271 -0.029397698 0.039469343 -0.1116858276
## RS6D -0.04656323 0.071344990 -0.282674982 0.034322064 0.1192637551
## RS7B -0.13664499 -0.108397584 0.272633413 0.224680179 -0.1356225698
## RS7C -0.36398426 -0.200868078 -0.047441111 0.181188310 0.0559780216
## RS8B 0.03695795 0.064064360 0.197500072 0.179495826 -0.1395988621
## RS8C 0.03472902 0.045298235 0.277202598 0.222178199 -0.1194076648
## RS9A 0.08963872 0.288165080 -0.078615357 -0.010370655 0.2943785816
## RS9C -0.02230774 -0.074912507 0.057974312 0.118356057 -0.0539003158
## MEpink MEpurple MEred MEsalmon MEtan
## RF13B -0.0647255585 0.0151669067 -0.0412048930 -0.029662215 0.004223004
## RF13D 0.0567711235 0.0045436211 0.0182869908 0.049450663 -0.066235131
## RF14B -0.0227504438 0.2937518176 0.0275893621 0.103027496 0.244085251
## RF14C -0.0264419838 0.2502514788 0.0597864506 0.001006515 0.166964022
## RF15B 0.2977573509 -0.1588172587 0.2684242914 0.307372343 -0.248831159
## RF15D 0.2675751469 -0.1318729337 0.2400972161 0.345022730 -0.142520988
## RF17B 0.0201348420 0.2433856181 0.0331804376 0.081061125 0.135198540
## RF17D -0.0450480509 0.3059433164 0.0755925194 -0.096919085 0.250243406
## RF18B 0.1219523070 0.2395830981 0.1781091787 -0.020943218 0.091557210
## RF18D 0.0163742929 0.2744014056 0.0723252141 -0.127858562 0.059736901
## RF19B -0.1720667112 0.2498582943 -0.0512888820 -0.108542763 0.116936546
## RF19C 0.1635415194 0.2285495607 0.1172686495 0.229733153 0.163346305
## RF20B 0.0648717076 0.2523383762 0.1687415292 -0.114089695 0.103224721
## RF20C -0.0006323075 0.2815220169 0.1134350350 0.015289757 0.220902418
## RF22B 0.0149928171 0.0064553830 0.0784595688 0.066325514 0.018193497
## RF22C 0.0863902854 0.0148348145 0.2755034633 -0.231270487 -0.002092525
## RF23A 0.0596190520 -0.1409510451 0.1214340863 0.105376690 -0.210702052
## RF23C 0.3621427737 -0.1708847815 0.3540544053 0.288780993 -0.283396351
## RF24B 0.0647018282 -0.1184864108 0.0693294328 0.104284706 -0.251013442
## RF24D 0.0353817295 -0.1166670846 0.1366884407 0.031922256 -0.114720296
## RF25A -0.0562018151 0.0289125326 -0.0087538500 0.068680685 0.067021221
## RF25C -0.1662642378 0.0155097767 -0.1071282788 -0.140434773 0.027334367
## RS11B 0.1710590129 -0.0851167145 0.0177296440 0.071706946 -0.079916232
## RS11D -0.0123017550 -0.0628538011 -0.0062733166 -0.076054829 -0.015876831
## RS12A -0.4251690594 -0.0525968459 -0.2413145580 -0.199216173 0.121356930
## RS12C -0.2233364331 -0.0952775719 -0.1958107921 0.158852236 0.209851744
## RS13A -0.0548748042 0.0002636796 -0.1072797695 -0.183481687 -0.095731690
## RS13C 0.1051871651 -0.0270560269 0.0968598729 -0.099625983 0.009833324
## RS14B -0.0590788210 -0.1035488583 -0.1273157643 0.042699790 0.078095475
## RS14C 0.0352559251 -0.1717060997 -0.1259982955 0.097974513 0.114293125
## RS15B 0.0304411988 -0.0565690422 -0.0572671524 -0.067166091 0.116585393
## RS15D -0.2099670058 -0.0853186762 -0.3806235298 0.071515739 0.114460742
## RS1B -0.2663953476 -0.1157331594 -0.2840088869 -0.330538257 -0.246867986
## RS1C 0.0357629036 -0.1631654026 -0.0767202606 0.029871371 -0.237046609
## RS2B 0.0606121877 -0.0915451370 -0.0867856299 0.140848204 -0.128993836
## RS2C 0.0099759699 -0.1103592577 -0.0290874453 -0.012361104 -0.183766930
## RS3B 0.0387810590 -0.0995603565 0.0381724975 0.007501047 -0.204687975
## RS3D 0.0481439887 -0.0967228437 0.0557197753 -0.014781091 -0.124814400
## RS6A 0.0257448243 -0.0793110778 -0.0310425919 -0.051437010 -0.088565740
## RS6D 0.1257730641 -0.0832343533 -0.0373659007 0.129638467 -0.088335417
## RS7B -0.1065909184 -0.0385199251 -0.1463727209 -0.183028241 0.161101074
## RS7C 0.0259720871 -0.1106157481 -0.1274810919 0.158323826 0.151502090
## RS8B -0.2639940985 -0.0451161373 -0.1818321106 -0.202533304 0.102274828
## RS8C -0.2604526033 -0.0217115573 -0.1850355018 -0.267425678 0.045501922
## RS9A 0.0598996333 -0.0546996378 0.0188980123 -0.093487238 -0.062180990
## RS9C 0.0314761592 -0.0172539519 0.0003051485 -0.055409283 -0.017527481
## MEturquoise MEyellow
## RF13B -0.049987446 0.063928602
## RF13D -0.019137147 -0.027995593
## RF14B -0.093185651 0.087636966
## RF14C -0.147709679 0.028957143
## RF15B 0.035232756 -0.382862198
## RF15D 0.078786749 -0.140150323
## RF17B -0.103854472 0.071280185
## RF17D -0.215289310 0.133312986
## RF18B -0.192378713 0.007020303
## RF18D -0.265041143 0.169001768
## RF19B -0.186757003 0.128738799
## RF19C 0.040882662 -0.155862228
## RF20B -0.249292334 0.114648118
## RF20C -0.150262215 0.146705924
## RF22B 0.091541606 0.028452390
## RF22C -0.186848856 0.090375080
## RF23A -0.089565674 0.014122502
## RF23C 0.013769414 -0.153623198
## RF24B -0.105806483 -0.039796849
## RF24D -0.105876790 -0.006773177
## RF25A 0.040721698 0.132236412
## RF25C -0.081198520 0.122990690
## RS11B 0.130518859 -0.158996790
## RS11D -0.014869673 0.054816786
## RS12A -0.007885866 0.293024916
## RS12C 0.308619894 -0.086322820
## RS13A -0.147909312 0.216107093
## RS13C -0.027925323 0.115750928
## RS14B 0.165844629 0.057513137
## RS14C 0.326288209 0.001611669
## RS15B 0.130767735 -0.019966361
## RS15D 0.307123864 -0.172427392
## RS1B -0.101737577 0.104667302
## RS1C 0.143128087 -0.183411878
## RS2B 0.100982719 -0.083769236
## RS2C 0.024999718 -0.048380219
## RS3B 0.024856931 -0.037476547
## RS3D 0.065598226 -0.130456228
## RS6A 0.076021058 -0.120889165
## RS6D 0.237997729 -0.426863340
## RS7B 0.010898755 0.148428399
## RS7C 0.300330368 -0.232336901
## RS8B -0.059359541 0.115675577
## RS8C -0.069784904 0.204699604
## RS9A -0.029634459 -0.091637783
## RS9C 0.046386423 0.048294943library(dplyr)
MEs <- MEs %>%
select_if(~ !any(is.na(.)))
# Calculate dissimilarity of module eigengenes
MEDiss = 1-cor(MEs, use = "pairwise.complete.obs");
# Cluster module eigengenes
METree = hclust(as.dist(MEDiss), method = "average");
# Plot the result
sizeGrWindow(7, 6)
png("../../output/WGCNA/ClusterEigengenes.png", 1000, 1000, pointsize=20)
plot(METree, main = "Clustering of module eigengenes - SIGNED (dissimilarity calc = MEDiss = 1-cor(MEs))",
xlab = "", sub = "")
MEDissThres = 0.15
abline(h=MEDissThres, col = "red")
dev.off()## pdf
## 2Step 7: Specify the cut line for the dendrogram (module) - Calc MODULE EIGENGENES (mergeMEs)
We had 17 modules before merging, and 15 modules after merging
MEDissThres = 0.15 # **Merge modules with >85% eigengene similarity.** Most studies use somewhere between 80-90% similarity. I will use 85% similarity as my merging threshold.
# Plot the cut line into the dendrogram
#abline(h=MEDissThres, col = "red")
# Call an automatic merging function
# merge = mergeCloseModules(dds.d0_vst, dynamicColors, cutHeight = MEDissThres, verbose = 3)
merge = mergeCloseModules(vst2, dynamicColors_sign, cutHeight = MEDissThres, verbose = 3)## mergeCloseModules: Merging modules whose distance is less than 0.15
## multiSetMEs: Calculating module MEs.
## Working on set 1 ...
## moduleEigengenes: Calculating 17 module eigengenes in given set.
## multiSetMEs: Calculating module MEs.
## Working on set 1 ...
## moduleEigengenes: Calculating 15 module eigengenes in given set.
## Calculating new MEs...
## multiSetMEs: Calculating module MEs.
## Working on set 1 ...
## moduleEigengenes: Calculating 15 module eigengenes in given set.# The merged module colors
mergedColors = merge$colors;
# Eigengenes of the new merged modules:
mergedMEs = merge$newMEs;
library(dplyr)
mergedMEs <- mergedMEs %>%
select_if(~ !any(is.na(.)))
# Cluster module eigengenes
MEDiss2 = 1-cor(mergedMEs,use = 'pairwise.complete.obs');
MEDiss2## MEcyan MEblack MEsalmon MEpink MEred MEgreen
## MEcyan 0.0000000 0.8363810 9.456342e-01 1.0357865 0.9335928 0.8183041
## MEblack 0.8363810 0.0000000 5.863287e-01 0.6853443 0.6615370 0.9676272
## MEsalmon 0.9456342 0.5863287 3.330669e-16 0.3469933 0.5950190 0.9701274
## MEpink 1.0357865 0.6853443 3.469933e-01 0.0000000 0.1847992 0.7966365
## MEred 0.9335928 0.6615370 5.950190e-01 0.1847992 0.0000000 0.5869122
## MEgreen 0.8183041 0.9676272 9.701274e-01 0.7966365 0.5869122 0.0000000
## MEmidnightblue 0.8698516 0.8018912 5.817853e-01 0.6693067 0.7309443 0.2592991
## MEblue 0.9280028 0.7794428 1.490247e+00 1.0050357 0.6357651 0.5714618
## MEbrown 0.8001074 0.9517543 7.364483e-01 0.5916456 0.3078619 0.4095681
## MEpurple 0.9256202 1.5772808 1.197132e+00 1.0905539 0.7866850 0.5726397
## MEturquoise 1.1349435 0.9426693 5.237657e-01 0.9541162 1.4074916 1.5287035
## MEtan 0.9540523 1.6478878 1.185038e+00 1.4159836 1.2577153 1.0651031
## MEmagenta 1.0943541 1.4972583 1.596656e+00 1.7187830 1.8343444 1.5388426
## MElightcyan 0.9911347 1.2290197 1.671253e+00 1.7962628 1.6277834 1.4917147
## MEyellow 0.8784473 1.2609555 1.706417e+00 1.5933035 1.2003612 1.0077827
## MEmidnightblue MEblue MEbrown MEpurple MEturquoise
## MEcyan 0.8698516 0.9280028 0.8001074 0.9256202 1.1349435
## MEblack 0.8018912 0.7794428 0.9517543 1.5772808 0.9426693
## MEsalmon 0.5817853 1.4902471 0.7364483 1.1971317 0.5237657
## MEpink 0.6693067 1.0050357 0.5916456 1.0905539 0.9541162
## MEred 0.7309443 0.6357651 0.3078619 0.7866850 1.4074916
## MEgreen 0.2592991 0.5714618 0.4095681 0.5726397 1.5287035
## MEmidnightblue 0.0000000 1.0911109 0.7673563 1.0220637 0.9399263
## MEblue 1.0911109 0.0000000 0.6094019 0.7720365 1.8569108
## MEbrown 0.7673563 0.6094019 0.0000000 0.2973886 1.6531319
## MEpurple 1.0220637 0.7720365 0.2973886 0.0000000 1.6328755
## MEturquoise 0.9399263 1.8569108 1.6531319 1.6328755 0.0000000
## MEtan 1.1727525 1.3358764 0.8763249 0.3437970 1.0458386
## MEmagenta 1.4884948 1.2900461 1.7130179 1.1285645 0.6777335
## MElightcyan 1.6092707 0.9317914 1.4922356 1.0135890 1.0509252
## MEyellow 1.4254228 0.5307596 0.9419223 0.5606714 1.6026729
## MEtan MEmagenta MElightcyan MEyellow
## MEcyan 0.9540523 1.0943541 9.911347e-01 0.8784473
## MEblack 1.6478878 1.4972583 1.229020e+00 1.2609555
## MEsalmon 1.1850376 1.5966556 1.671253e+00 1.7064168
## MEpink 1.4159836 1.7187830 1.796263e+00 1.5933035
## MEred 1.2577153 1.8343444 1.627783e+00 1.2003612
## MEgreen 1.0651031 1.5388426 1.491715e+00 1.0077827
## MEmidnightblue 1.1727525 1.4884948 1.609271e+00 1.4254228
## MEblue 1.3358764 1.2900461 9.317914e-01 0.5307596
## MEbrown 0.8763249 1.7130179 1.492236e+00 0.9419223
## MEpurple 0.3437970 1.1285645 1.013589e+00 0.5606714
## MEturquoise 1.0458386 0.6777335 1.050925e+00 1.6026729
## MEtan 0.0000000 0.6018516 6.027920e-01 0.5881998
## MEmagenta 0.6018516 0.0000000 2.388853e-01 0.6635948
## MElightcyan 0.6027920 0.2388853 1.110223e-16 0.2320545
## MEyellow 0.5881998 0.6635948 2.320545e-01 0.0000000METree2 = hclust(as.dist(MEDiss2), method = "average");
# Plot the result
plot(METree2, main = "Clustering of module eigengenes - SIGNED (dissimilarity calc = MEDiss = 1-cor(MEs))",
xlab = "", sub = "")
sizeGrWindow(7, 6)
png("../../output/WGCNA/ClusterEigengenes_merged.png", 1000, 1000, pointsize=20)
plot(METree2, main = "Clustering of module eigengenes - SIGNED (dissimilarity calc = MEDiss = 1-cor(MEs))",
xlab = "", sub = "")
dev.off()## pdf
## 2Step 8: Plot dendrogram with the cut line ‘MEDissThres’
plotDendroAndColors(geneTree_sign, cbind(dynamicColors_sign, mergedColors),
c("Dynamic Tree Cut", "Merged dynamic"),
dendroLabels = FALSE, hang = 0.03,
addGuide = TRUE, guideHang = 0.05)
sizeGrWindow(12, 9)
png("../../output/WGCNA/ClusterDendrogram_signed.png", 1000, 1000, pointsize=20)
plotDendroAndColors(geneTree_sign, cbind(dynamicColors_sign, mergedColors),
c("Dynamic Tree Cut", "Merged dynamic"),
dendroLabels = FALSE, hang = 0.03,
addGuide = TRUE, guideHang = 0.05)
dev.off()## pdf
## 2Step 9: Commit to mergedcolors as ‘MEs’ and ‘moduleColors’
# Rename to moduleColors
moduleColors = mergedColors
# Construct numerical labels corresponding to the colors
colorOrder = c("grey", standardColors(50));
moduleLabels = match(moduleColors, colorOrder)-1;
MEs = mergedMEs;
# Save module colors and labels for use in subsequent parts
save(MEs, moduleLabels, moduleColors, geneTree_sign, file = "../../output/WGCNA/networkConstruction-stepByStep.RData")
# write csv - save the module eigengenes
write.csv(MEs, file = "../../output/WGCNA/WGCNA_ModuleEigengenes.csv")
table(mergedColors)## mergedColors
## black blue brown cyan green lightcyan
## 396 1989 942 105 563 65
## magenta midnightblue pink purple red salmon
## 219 71 220 218 425 154
## tan turquoise yellow
## 170 2558 917Prepare for module trait associations - Eigengene calc - trait data as factors
#Prepare trait data. Data has to be numeric, so I will substitute time_points and type for numeric values. The "trait" we are considering here is habitat of origin. Make a dataframe that has a column for each origin name and a row for samples. Populate a 1 for samples that match each origin = flat and a 0 for samples not matching (e.g., origin = slope). This process changes origin from a categorical variable into a binary variable. This will allow for correlations between mean eigengenes and habitat of origin.
Samples = rownames(vst2);# start new variable 'd21.Samples' calling the row names of the gene data (sample as 'Mouse' in trait data)
TreatRows = match(Samples, coldata$Coral_ID); # match the names
Traits = coldata[TreatRows, -1]; # removes the row numbers
rownames(Traits) = coldata[TreatRows, 1]; # inserts the new traitRows - matches sample treatments
dim(Traits) # 46 Samples 2 columns (Origin and Treatment)## [1] 46 3all(rownames(Traits) == rownames(vst2)) # should be TRUE## [1] TRUEdim(Traits) # 46 2## [1] 46 3Origin
# ALL TRAITS
head(Traits)## Origin Treatment Colony
## RF13B Flat Variable RF13
## RF13D Flat Stable RF13
## RF14B Flat Variable RF14
## RF14C Flat Stable RF14
## RF15B Flat Variable RF15
## RF15D Flat Stable RF15# Origin groups
Traits.Origin <- Traits %>% dplyr::select('Origin')
Traits.Origin$Flat <- (Traits.Origin$Origin == "Flat")
Traits.Origin$Flat <- as.numeric(Traits.Origin$Flat)
Traits.Origin$Flat <- as.factor(Traits.Origin$Flat)
Traits.Origin$Slope <- (Traits.Origin$Origin == "Slope")
Traits.Origin$Slope <- as.numeric(Traits.Origin$Slope)
Traits.Origin$Slope <- as.factor(Traits.Origin$Slope)
Traits.Origin <- Traits.Origin[,c(2:3)]
head(Traits.Origin)## Flat Slope
## RF13B 1 0
## RF13D 1 0
## RF14B 1 0
## RF14C 1 0
## RF15B 1 0
## RF15D 1 0Cluster samples by Origin
# Primary treatment Only
traitColors_Primary = labels2colors(Traits.Origin); # Convert traits to a color representation: white means low, red means high, grey means missing entry
plotDendroAndColors(sampleTree, traitColors_Primary, # Plot the sample dendrogram and the colors underneath.
groupLabels = names(Traits.Origin),
main = "Sample dendrogram and trait heatmap")
png("../../output/WGCNA/ClusterTree_Origin.png", 1000, 1000, pointsize=20)
traitColors_Primary = labels2colors(Traits.Origin); # Convert traits to a color representation: white means low, red means high, grey means missing entry
plotDendroAndColors(sampleTree, traitColors_Primary, # Plot the sample dendrogram and the colors underneath.
groupLabels = names(Traits.Origin),
main = "Sample dendrogram and trait heatmap")
dev.off()## quartz_off_screen
## 2# save data
save(vst2, Traits, file = "../../output/WGCNA/dataInput.RData")
# write the vst transformed data
write.csv(vst2, "../../output/WGCNA/vstTransformed_WGCNAdata.csv") # write# identify modules that are significantly associated with the measured clinical traits.
# Since we already have a summary profile (eigengene) for each module, we simply correlate eigengenes with external traits and look for the most significant associations:
# Define numbers of genes and samples
nGenes = ncol(vst2); # 9266
nSamples = nrow(vst2); # 46
# # Recalculate MEs with color labels
MEs0 = moduleEigengenes(vst2, moduleColors)$eigengenes
MEs = orderMEs(MEs0) # reorders the columns (colors/modules) , doesn't change anything
# change character treatments to integers
# ALL TRAIT DATA
Traits$Origin <- as.factor(Traits$Origin)
Traits$Origin <- as.numeric(Traits$Origin)Module trait correlation
# ALL TRAIT DATA
dim(Traits) # 48 2## [1] 46 3dim(MEs) # 48 15## [1] 46 15# moduleTraitCor = cor(MEs, d0.Traits, use = "p");
# moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples);
nSamples = nrow(vst2)
# Origin
Traits.Origin$Flat <- as.numeric(Traits.Origin$Flat)
Traits.Origin$Slope <- as.numeric(Traits.Origin$Slope)
moduleTraitCor_Primary = cor(MEs, Traits.Origin, use = "p");
moduleTraitPvalue_Primary = corPvalueStudent(moduleTraitCor_Primary, nSamples);Heatmaps
sizeGrWindow(10,10)
# Will display correlations and their p-values
d0.PRIMARYTreatments.matrix <- paste(signif(moduleTraitCor_Primary, 3), "\n(",
signif(moduleTraitPvalue_Primary, 3), ")", sep = "")
#dim(textMatrix) == dim(moduleTraitCor_treatonly)
par(mar = c(8, 9.5, 5, 3));
# Display the correlation values within a heatmap plot
labeledHeatmap(Matrix = moduleTraitCor_Primary,
xLabels = names(Traits.Origin),
yLabels = names(MEs),
ySymbols = names(MEs),
colorLabels = TRUE,
colors = blueWhiteRed(50),
textMatrix = d0.PRIMARYTreatments.matrix,
setStdMargins = FALSE,
cex.text = 1,
zlim = c(-1,1),
main = paste("Module-trait relationships - Origin"))
png("../../output/WGCNA/Treatments_Primary_heatmap2.png", 500, 1000, pointsize=20)
labeledHeatmap(Matrix = moduleTraitCor_Primary,
xLabels = names(Traits.Origin),
yLabels = names(MEs),
ySymbols = names(MEs),
colorLabels = TRUE,
colors = blueWhiteRed(50),
textMatrix = d0.PRIMARYTreatments.matrix,
setStdMargins = FALSE,
cex.text = 1,
zlim = c(-1,1),
main = paste("Module-trait relationships - Origin"))
dev.off()## pdf
## 2Module eigengene - MEs boxplots by treatment group
library(reshape2) ##
## Attaching package: 'reshape2'## The following object is masked from 'package:tidyr':
##
## smithsMEs_table <- MEs # new table for plotting
MEs_table$Coral_ID <- rownames(MEs) # call rows as coolumn to merge with treatment data
MEsPlotting <- merge(coldata, MEs_table, by = 'Coral_ID') %>% dplyr::select(-c("Colony")) # merge
#MEsPlotting <- MEsPlotting[,-c(2)] # ommit the phys data to just plot the module colors
MEsPlotting_melt <- melt(MEsPlotting, id=c('Coral_ID', 'Origin', 'Treatment'))
#plot it
ggplot(MEsPlotting_melt, aes(x=Origin, y=value, fill = factor(Origin), shape=Origin)) +
geom_boxplot(aes(middle = mean(value)), position=position_dodge(0.8), outlier.size = 0, alpha = 0.5) +
stat_summary(fun = mean, color = "black", position = position_dodge(0.75),
geom = "point", shape = 19, size = 3,
show.legend = FALSE) +
ylab("ModuleEigengene") +
ylim(-0.5,0.5) +
scale_color_manual(values=c("#56B4E9","#D55E00")) +
geom_hline(yintercept=0, linetype='dotted', col = 'black', size = 1)+
theme_classic() +
theme(legend.position = "none") +
facet_wrap(~variable)## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_boxplot()`).## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_summary()`).
png("../../output/WGCNA/ME_Boxplot.png", 600, 1000, pointsize=20)
ggplot(MEsPlotting_melt, aes(x=Origin, y=value, fill = factor(Origin), shape=Origin)) +
geom_boxplot(aes(middle = mean(value)), position=position_dodge(0.8), outlier.size = 0, alpha = 0.5) +
stat_summary(fun = mean, color = "black", position = position_dodge(0.75),
geom = "point", shape = 19, size = 3,
show.legend = FALSE) +
ylab("ModuleEigengene") +
ylim(-0.5,0.5) +
scale_color_manual(values=c("#56B4E9","#D55E00")) +
geom_hline(yintercept=0, linetype='dotted', col = 'black', size = 1)+
theme_classic() +
theme(legend.position = "none") +
facet_wrap(~variable)## Warning: Removed 1 row containing non-finite outside the scale range (`stat_boxplot()`).
## Removed 1 row containing non-finite outside the scale range (`stat_summary()`).dev.off()## quartz_off_screen
## 2Module trait correlation
Load in phys for all traits
d <- read.csv("../../data/Heron pHi coral physiology and respirometry R_reduced.csv", strip.white=T)
d$Origin <- as.factor(d$Origin)
d$Treatment <- as.factor(d$Treatment)
d$Origin <- as.numeric(d$Origin)
d$Treatment <- as.numeric(d$Treatment)
d$Colony <- as.factor(d$Colony)
d$Coral.ID <- as.factor(d$Coral.ID)
d <- plyr::rename(d, c("Coral.ID"="Coral_ID"))
d <- d[-c(7, 8),]#here we remove the two samples that are behaving weirdly
d <- subset(d, select = -c(Coral_ID))Correlations of traits with eigengenes
moduleTraitCor_both = cor(MEs, d, use = "p")## Warning in storage.mode(y) <- "double": NAs introduced by coercionmoduleTraitPvalue_both = corPvalueStudent(moduleTraitCor_both, nSamples)
Colors=sub("ME","",names(MEs))
moduleTraitTree = hclust(dist(t(moduleTraitCor_both)), method = "average")
plot(moduleTraitTree, main = "Clustering based on module-trait correlation", sub="", xlab="", cex.lab = 1.5, cex.axis = 1.5, cex.main = 2)
### Heatmaps
#sizeGrWindow(20,20)
# Will display correlations and their p-values
both.matrix <- paste(signif(moduleTraitCor_both, 3), "\n(",
signif(moduleTraitPvalue_both, 3), ")", sep = "")
#dim(textMatrix) == dim(moduleTraitCor_treatonly)
sizeGrWindow(20,20)
labeledHeatmap(Matrix = moduleTraitCor_both,
xLabels = names(d),
yLabels = names(MEs),
ySymbols = names(MEs),
colorLabels = TRUE,
colors = blueWhiteRed(50),
textMatrix = both.matrix,
setStdMargins = FALSE,
cex.text = 1,
zlim = c(-1,1),
main = paste("Module-trait relationships - Origin and Treatment"))
par(mar = c(10, 9.5, 5, 3));
# Display the correlation values within a heatmap plot
png("../../output/WGCNA/Both_with phys and pHi_heatmap_new.png", 2000, 2000, pointsize=20)
labeledHeatmap(Matrix = moduleTraitCor_both,
xLabels = names(d),
yLabels = names(MEs),
ySymbols = names(MEs),
colorLabels = TRUE,
colors = blueWhiteRed(50),
textMatrix = both.matrix,
setStdMargins = FALSE,
cex.text = 1,
zlim = c(-1,1),
main = paste("Module-trait relationships - Origin and Treatment"))
dev.off()## pdf
## 2Complex heatmap
library(ComplexHeatmap)## Loading required package: grid## ========================================
## ComplexHeatmap version 2.16.0
## Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
## Github page: https://github.com/jokergoo/ComplexHeatmap
## Documentation: http://jokergoo.github.io/ComplexHeatmap-reference
##
## If you use it in published research, please cite either one:
## - Gu, Z. Complex Heatmap Visualization. iMeta 2022.
## - Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional
## genomic data. Bioinformatics 2016.
##
##
## The new InteractiveComplexHeatmap package can directly export static
## complex heatmaps into an interactive Shiny app with zero effort. Have a try!
##
## This message can be suppressed by:
## suppressPackageStartupMessages(library(ComplexHeatmap))
## ========================================
## ! pheatmap() has been masked by ComplexHeatmap::pheatmap(). Most of the arguments
## in the original pheatmap() are identically supported in the new function. You
## can still use the original function by explicitly calling pheatmap::pheatmap().##
## Attaching package: 'ComplexHeatmap'## The following object is masked from 'package:pheatmap':
##
## pheatmap## The following object is masked from 'package:genefilter':
##
## dist2#bold sig p-values
#dendrogram with WGCNA MEtree cut-off
#colored y-axis
#Create list of pvalues for eigengene correlation with specific life stages
heatmappval <- signif(moduleTraitPvalue_both, 1)
#Make list of heatmap row colors
htmap.colors <- rownames(moduleTraitCor_both)
htmap.colors <- gsub("ME", "", htmap.colors)
moduleTraitCor_both <- as.matrix(moduleTraitCor_both)
row_dend = as.dendrogram(hclust(dist(moduleTraitCor_both), method = "average"))
row_dend = rotate(row_dend,c(1:3,5,4,8:10, 12, 11, 14:15 ,13 ,6:7))
column_dend = as.dendrogram(hclust(dist(t(moduleTraitCor_both)), method = "average"))
column_dend = rotate(column_dend,c(1:2,12:17,7:11,3:6))
pdf("../../output/WGCNA/Both_with phys and pHi_heatmap_new_row_clust.pdf", height = 40, width = 40)
ht=Heatmap(moduleTraitCor_both, name = "Eigengene", column_title = "Module-Trait Eigengene Correlation",
col = blueWhiteRed(50),
row_names_side = "left",
cluster_rows = row_dend,
row_dend_width = unit(2, "in"),
width = unit(20, "in"), height = unit(20, "in"),
cluster_columns = column_dend,
column_dend_height = unit(1.25, "in"),
row_gap = unit(2.5, "mm"), border = TRUE,
cell_fun = function(j, i, x, y, w, h, col) {
if(heatmappval[i, j] < 0.05) {
grid.text(sprintf("%s", heatmappval[i, j]), x, y, gp = gpar(fontsize = 25, fontface = "bold"))
}
else {
#grid.text(sprintf("%s", heatmappval[i, j]), x, y, gp = gpar(fontsize = 25, fontface = "plain"))
}},
column_names_gp = gpar(fontsize = 30), row_names_gp = gpar(fontsize = 25, alpha = 0.75, border = TRUE, fill = htmap.colors))
draw(ht)
dev.off()## quartz_off_screen
## 2Module eigengene - MEs boxplots by treatment group
library(reshape2)
MEs_table <- MEs # new table for plotting
MEs_table$Coral_ID <- rownames(MEs) # call rows as coolumn to merge with treatment data
MEsPlotting <- merge(coldata, MEs_table, by = 'Coral_ID') %>% dplyr::select(-c("Colony")) # merge
#MEsPlotting <- MEsPlotting[,-c(2)] # ommit the phys data to just plot the module colors
MEsPlotting_melt <- melt(MEsPlotting, id=c('Coral_ID', 'Origin', 'Treatment'))
#plot it
png("../../output/WGCNA/ME_Boxplot_Treatment.png", 600, 1000, pointsize=20)
ME_Boxplot_Treatment = ggplot(MEsPlotting_melt, aes(x=Treatment, y=value, fill = factor(Treatment), shape=Treatment)) +
geom_boxplot(aes(middle = mean(value)), position=position_dodge(0.8), outlier.size = 0, alpha = 0.5) +
stat_summary(fun.y = mean, color = "black", position = position_dodge(0.75),
geom = "point", shape = 19, size = 3,
show.legend = FALSE) +
ylab("ModuleEigengene") +
ylim(-0.5,0.5) +
scale_color_manual(values=c("#56B4E9","#D55E00")) +
geom_hline(yintercept=0, linetype='dotted', col = 'black', size = 1)+
theme_classic() +
theme(legend.position = "none") +
facet_wrap(~variable)## Warning: The `fun.y` argument of `stat_summary()` is deprecated as of ggplot2 3.3.0.
## ℹ Please use the `fun` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.ME_Boxplot_Treatment## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_boxplot()`).## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_summary()`).dev.off()## quartz_off_screen
## 2ME_Boxplot_Treatment## Warning: Removed 1 row containing non-finite outside the scale range (`stat_boxplot()`).
## Removed 1 row containing non-finite outside the scale range (`stat_summary()`).
Plots based off of Strader & Quigley 2022: https://www.nature.com/articles/s41467-022-32217-z
# View module eigengene data and make dataframe for Strader plots.
head(MEs)## MEcyan MEblack MEsalmon MEpink MEred MEgreen
## RF13B 0.023045712 0.03573889 -0.029662215 -0.06472556 -0.04120489 0.06229807
## RF13D 0.010859993 0.01671523 0.049450663 0.05677112 0.01828699 0.04001668
## RF14B 0.008973815 -0.14656653 0.103027496 -0.02275044 0.02758936 0.09009377
## RF14C 0.055555248 -0.14441784 0.001006515 -0.02644198 0.05978645 0.12396450
## RF15B 0.066984890 0.36136333 0.307372343 0.29775735 0.26842429 0.07146349
## RF15D 0.061369555 0.34069144 0.345022730 0.26757515 0.24009722 0.12177161
## MEmidnightblue MEblue MEbrown MEpurple MEturquoise
## RF13B -0.024549173 0.005857733 0.08621400 0.015166907 -0.04998745
## RF13D -0.018376306 0.044219420 0.07779582 0.004543621 -0.01913715
## RF14B -0.022880675 -0.083031859 0.24394568 0.293751818 -0.09318565
## RF14C 0.006327375 0.035898347 0.23139138 0.250251479 -0.14770968
## RF15B 0.151329333 0.033563498 0.19519443 -0.158817259 0.03523276
## RF15D 0.259392975 -0.059028765 0.13060026 -0.131872934 0.07878675
## MEtan MEmagenta MElightcyan MEyellow
## RF13B 0.004223004 0.040533178 0.004539119 0.06392860
## RF13D -0.066235131 0.004749849 -0.048921622 -0.02799559
## RF14B 0.244085251 -0.055936540 0.007711938 0.08763697
## RF14C 0.166964022 -0.015987095 -0.069472828 0.02895714
## RF15B -0.248831159 -0.362117043 -0.344930466 -0.38286220
## RF15D -0.142520988 -0.295093407 -0.167570877 -0.14015032names(MEs)## [1] "MEcyan" "MEblack" "MEsalmon" "MEpink"
## [5] "MEred" "MEgreen" "MEmidnightblue" "MEblue"
## [9] "MEbrown" "MEpurple" "MEturquoise" "MEtan"
## [13] "MEmagenta" "MElightcyan" "MEyellow"Strader_MEs <- MEs
Strader_MEs$Treatment <- coldata_outrm$Treatment
Strader_MEs$Origin <- coldata_outrm$Origin
Strader_MEs$Colony <- d$Colony
Strader_MEs$Coral_ID <- rownames(Strader_MEs)
head(Strader_MEs)## MEcyan MEblack MEsalmon MEpink MEred MEgreen
## RF13B 0.023045712 0.03573889 -0.029662215 -0.06472556 -0.04120489 0.06229807
## RF13D 0.010859993 0.01671523 0.049450663 0.05677112 0.01828699 0.04001668
## RF14B 0.008973815 -0.14656653 0.103027496 -0.02275044 0.02758936 0.09009377
## RF14C 0.055555248 -0.14441784 0.001006515 -0.02644198 0.05978645 0.12396450
## RF15B 0.066984890 0.36136333 0.307372343 0.29775735 0.26842429 0.07146349
## RF15D 0.061369555 0.34069144 0.345022730 0.26757515 0.24009722 0.12177161
## MEmidnightblue MEblue MEbrown MEpurple MEturquoise
## RF13B -0.024549173 0.005857733 0.08621400 0.015166907 -0.04998745
## RF13D -0.018376306 0.044219420 0.07779582 0.004543621 -0.01913715
## RF14B -0.022880675 -0.083031859 0.24394568 0.293751818 -0.09318565
## RF14C 0.006327375 0.035898347 0.23139138 0.250251479 -0.14770968
## RF15B 0.151329333 0.033563498 0.19519443 -0.158817259 0.03523276
## RF15D 0.259392975 -0.059028765 0.13060026 -0.131872934 0.07878675
## MEtan MEmagenta MElightcyan MEyellow Treatment Origin
## RF13B 0.004223004 0.040533178 0.004539119 0.06392860 Variable Flat
## RF13D -0.066235131 0.004749849 -0.048921622 -0.02799559 Stable Flat
## RF14B 0.244085251 -0.055936540 0.007711938 0.08763697 Variable Flat
## RF14C 0.166964022 -0.015987095 -0.069472828 0.02895714 Stable Flat
## RF15B -0.248831159 -0.362117043 -0.344930466 -0.38286220 Variable Flat
## RF15D -0.142520988 -0.295093407 -0.167570877 -0.14015032 Stable Flat
## Colony Coral_ID
## RF13B 13 RF13B
## RF13D 13 RF13D
## RF14B 14 RF14B
## RF14C 14 RF14C
## RF15B 15 RF15B
## RF15D 15 RF15DStrader_MEs <- Strader_MEs%>%
droplevels() #drop unused levelplot_MEs <- Strader_MEs%>%
gather(., key="Module", value="Mean", 1:(ncol(Strader_MEs)-4))
head(plot_MEs)## Treatment Origin Colony Coral_ID Module Mean
## 1 Variable Flat 13 RF13B MEcyan 0.023045712
## 2 Stable Flat 13 RF13D MEcyan 0.010859993
## 3 Variable Flat 14 RF14B MEcyan 0.008973815
## 4 Stable Flat 14 RF14C MEcyan 0.055555248
## 5 Variable Flat 15 RF15B MEcyan 0.066984890
## 6 Stable Flat 15 RF15D MEcyan 0.061369555plot_MEs_mean <- summarySE(plot_MEs, measurevar='Mean', groupvars=c('Treatment','Origin', "Module"), na.rm=TRUE, conf.interval = 0.95)
plot_MEs_mean## Treatment Origin Module N Mean sd se
## 1 Stable Slope MEblack 12 -0.039754264 0.06202534 0.017905174
## 2 Stable Slope MEblue 12 -0.110056761 0.20407118 0.058910276
## 3 Stable Slope MEbrown 12 -0.122307021 0.06850331 0.019775203
## 4 Stable Slope MEcyan 12 0.013369322 0.04588326 0.013245356
## 5 Stable Slope MEgreen 12 -0.057759479 0.18209748 0.052567016
## 6 Stable Slope MElightcyan 12 0.008978748 0.14457269 0.041734541
## 7 Stable Slope MEmagenta 12 0.094525394 0.07778958 0.022455917
## 8 Stable Slope MEmidnightblue 12 0.022929424 0.17937023 0.051779724
## 9 Stable Slope MEpink 12 -0.024042545 0.13090715 0.037789639
## 10 Stable Slope MEpurple 12 -0.087106274 0.05005723 0.014450276
## 11 Stable Slope MEred 12 -0.084292611 0.12967533 0.037434044
## 12 Stable Slope MEsalmon 12 0.010043182 0.12510710 0.036115309
## 13 Stable Slope MEtan 12 -0.001827060 0.13790126 0.039808664
## 14 Stable Slope MEturquoise 12 0.137324385 0.15076081 0.043520897
## 15 Stable Slope MEyellow 12 -0.071252071 0.17217664 0.049703116
## 16 Stable Flat MEblack 11 0.028212561 0.20306332 0.061225894
## 17 Stable Flat MEblue 11 0.051722038 0.09691348 0.029220513
## 18 Stable Flat MEbrown 11 0.140572527 0.08941799 0.026960539
## 19 Stable Flat MEcyan 11 0.030429746 0.02490805 0.007510061
## 20 Stable Flat MEgreen 11 0.092491604 0.09756417 0.029416704
## 21 Stable Flat MElightcyan 11 -0.063292592 0.11152780 0.033626897
## 22 Stable Flat MEmagenta 11 -0.101329792 0.12740965 0.038415455
## 23 Stable Flat MEmidnightblue 11 0.051785345 0.11968160 0.036085359
## 24 Stable Flat MEpink 11 0.068162754 0.14864751 0.044818912
## 25 Stable Flat MEpurple 11 0.086921017 0.18455066 0.055644117
## 26 Stable Flat MEred 11 0.123264555 0.12782853 0.038541753
## 27 Stable Flat MEsalmon 11 0.033156651 0.18594831 0.056065525
## 28 Stable Flat MEtan 11 0.025414739 0.16774005 0.050575528
## 29 Stable Flat MEturquoise 11 -0.094356803 0.11144867 0.033603039
## 30 Stable Flat MEyellow 11 0.018812643 0.12508447 0.037714387
## 31 Variable Slope MEblack 12 -0.021871855 0.05108377 0.014746614
## 32 Variable Slope MEblue 12 0.023848447 0.10688538 0.030855151
## 33 Variable Slope MEbrown 12 -0.125277907 0.04818037 0.013908476
## 34 Variable Slope MEcyan 12 -0.067307852 0.28331241 0.081785248
## 35 Variable Slope MEgreen 12 -0.076561117 0.15351021 0.044314581
## 36 Variable Slope MElightcyan 12 0.097493355 0.16605295 0.047935357
## 37 Variable Slope MEmagenta 12 0.098641298 0.10807566 0.031198757
## 38 Variable Slope MEmidnightblue 12 -0.090179594 0.16882156 0.048734588
## 39 Variable Slope MEpink 12 -0.065797094 0.17240455 0.049768905
## 40 Variable Slope MEpurple 12 -0.068504434 0.03317717 0.009577425
## 41 Variable Slope MEred 12 -0.099034919 0.10326861 0.029811079
## 42 Variable Slope MEsalmon 12 -0.087344334 0.13762115 0.039727804
## 43 Variable Slope MEtan 12 -0.027294229 0.13771990 0.039756309
## 44 Variable Slope MEturquoise 12 0.024446994 0.09902844 0.028587048
## 45 Variable Slope MEyellow 12 0.035223379 0.14284911 0.041236987
## 46 Variable Flat MEblack 11 0.039015933 0.21596847 0.065116943
## 47 Variable Flat MEblue 11 0.042323395 0.11000329 0.033167240
## 48 Variable Flat MEbrown 11 0.129520122 0.08552768 0.025787566
## 49 Variable Flat MEcyan 11 0.028412288 0.03483993 0.010504633
## 50 Variable Flat MEgreen 11 0.054039955 0.06651343 0.020054555
## 51 Variable Flat MElightcyan 11 -0.052858793 0.12351732 0.037241874
## 52 Variable Flat MEmagenta 11 -0.109397509 0.13192715 0.039777533
## 53 Variable Flat MEmidnightblue 11 0.021578478 0.07129269 0.021495555
## 54 Variable Flat MEpink 11 0.029844125 0.11991092 0.036154503
## 55 Variable Flat MEpurple 11 0.082836119 0.17714043 0.053409849
## 56 Variable Flat MEred 11 0.076729115 0.09973119 0.030070085
## 57 Variable Flat MEsalmon 11 0.051171879 0.11835171 0.035684383
## 58 Variable Flat MEtan 11 0.006353940 0.16860756 0.050837094
## 59 Variable Flat MEturquoise 11 -0.082121065 0.10597943 0.031954002
## 60 Variable Flat MEyellow 11 0.020491385 0.14445032 0.043553411
## ci
## 1 0.03940902
## 2 0.12966064
## 3 0.04352493
## 4 0.02915283
## 5 0.11569922
## 6 0.09185711
## 7 0.04942514
## 8 0.11396640
## 9 0.08317443
## 10 0.03180484
## 11 0.08239177
## 12 0.07948926
## 13 0.08761828
## 14 0.09578885
## 15 0.10939582
## 16 0.13641979
## 17 0.06510736
## 18 0.06007182
## 19 0.01673346
## 20 0.06554450
## 21 0.07492540
## 22 0.08559497
## 23 0.08040319
## 24 0.09986276
## 25 0.12398282
## 26 0.08587638
## 27 0.12492177
## 28 0.11268930
## 29 0.07487224
## 30 0.08403289
## 31 0.03245708
## 32 0.06791173
## 33 0.03061235
## 34 0.18000812
## 35 0.09753574
## 36 0.10550501
## 37 0.06866800
## 38 0.10726410
## 39 0.10954062
## 40 0.02107977
## 41 0.06561374
## 42 0.08744031
## 43 0.08750305
## 44 0.06291967
## 45 0.09076200
## 46 0.14508959
## 47 0.07390122
## 48 0.05745828
## 49 0.02340578
## 50 0.04468433
## 51 0.08298007
## 52 0.08862987
## 53 0.04789508
## 54 0.08055725
## 55 0.11900456
## 56 0.06700033
## 57 0.07950976
## 58 0.11327210
## 59 0.07119795
## 60 0.09704305By Origin, Plot mean module eigengene for each module.
png("../../output/WGCNA/plot_MEs_fig_Origin.png", 600, 1000, pointsize=20)
plot_MEs_fig_Origin <- ggplot(plot_MEs, aes(y=Mean, x=Origin, color=Origin, fill=Origin))+
geom_point(alpha=0.8,position=position_jitterdodge(0.05))+
geom_boxplot(alpha=0.4, color="black", outlier.shape = NA)+
#geom_smooth(method=gam)+
geom_hline(yintercept = 0, linetype="dashed", color = 'black', size=0.7, show.legend = TRUE)+
facet_wrap(~Module)+
#scale_color_brewer(palette = "RdBu", direction=-1)+
#scale_fill_brewer(palette = "RdBu", direction=-1)+
#scale_x_continuous(expression(Temperature~(degree~C)), limits=c(23,38),breaks=c(24,26,28,30,32,34,36,38),expand = c(0.005, 0.005))+
scale_y_continuous(expression(Mean~module~eigenegene), limits=c(-0.5, 0.5), expand = c(0.005, 0.005))+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5, size=12),
axis.text.y=element_text(vjust=0.5, size=12),
axis.title.x=element_text(size=12),
axis.title.y=element_text(size=12),
legend.text = element_text(vjust=0.5, size=12),
panel.background= element_rect(fill=NA, color='black'),
#legend.position= c(0.89,0.89),
strip.text = element_text(vjust=0.5, size=12))
plot_MEs_fig_Origin
dev.off()## quartz_off_screen
## 2By Treatment, Plot mean module eigengene for each module.
png("../../output/WGCNA/plot_MEs_fig_Treatment.png", 600, 1000, pointsize=20)
plot_MEs_fig_Treatment <- ggplot(plot_MEs, aes(y=Mean, x=Treatment, color=Treatment, fill=Treatment))+
geom_point(alpha=0.8,position=position_jitterdodge(0.05))+
geom_boxplot(alpha=0.4, color="black", outlier.shape = NA)+
#geom_smooth(method=gam)+
geom_hline(yintercept = 0, linetype="dashed", color = 'black', size=0.7, show.legend = TRUE)+
facet_wrap(~Module)+
scale_color_manual("Treatment", values=c("Stable"= 'lightcyan',"Variable"='orange'))+
scale_fill_manual("Treatment", values=c("Stable"= 'lightcyan',"Variable"='orange'))+
#scale_x_continuous(expression(Temperature~(degree~C)), limits=c(23,38),breaks=c(24,26,28,30,32,34,36,38),expand = c(0.005, 0.005))+
scale_y_continuous(expression(Mean~module~eigenegene), limits=c(-0.5, 0.5), expand = c(0.005, 0.005))+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5, size=12),
axis.text.y=element_text(vjust=0.5, size=12),
axis.title.x=element_text(size=12),
axis.title.y=element_text(size=12),
legend.text = element_text(vjust=0.5, size=12),
panel.background= element_rect(fill=NA, color='black'),
#legend.position= c(0.89,0.89),
strip.text = element_text(vjust=0.5, size=12))
plot_MEs_fig_Treatment
dev.off()## quartz_off_screen
## 2By Treatment and Origin, Plot mean module eigengene for each module.
png("../../output/WGCNA/plot_MEs_fig_Trt_Orig.png", 600, 1000, pointsize=20)
plot_MEs_fig_Trt_Orig <- ggplot(plot_MEs, aes(y=Mean, x=Treatment, color=Origin, fill=Origin))+
geom_point(alpha=0.8,position=position_jitterdodge(0.05))+
geom_boxplot(alpha=0.4, color="black", outlier.shape = NA)+
#geom_smooth(method=gam)+
geom_hline(yintercept = 0, linetype="dashed", color = 'black', size=0.7, show.legend = TRUE)+
facet_wrap(~Module)+
scale_fill_manual("Origin", values=c("Slope"='paleturquoise3', "Flat"= "indianred"))+
scale_color_manual("Origin", values=c("Slope"='paleturquoise3', "Flat"= "indianred"))+
#scale_x_continuous(expression(Temperature~(degree~C)), limits=c(23,38),breaks=c(24,26,28,30,32,34,36,38),expand = c(0.005, 0.005))+
scale_y_continuous(expression(Mean~module~eigenegene), limits=c(-0.5, 0.5), expand = c(0.005, 0.005))+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5, size=12),
axis.text.y=element_text(vjust=0.5, size=12),
axis.title.x=element_text(size=12),
axis.title.y=element_text(size=12),
legend.text = element_text(vjust=0.5, size=12),
panel.background= element_rect(fill=NA, color='black'),
#legend.position= c(0.89,0.89),
strip.text = element_text(vjust=0.5, size=12))
plot_MEs_fig_Trt_Orig
dev.off()## quartz_off_screen
## 2geneModuleMembership - geneTraitSignificance - GSPvalue
# Gene relationship to trait and important modules:
# Gene Significance and Module Membership
# We quantify associations of individual genes with our trait of interest (TAOC)
# names (colors) of the modules
modNames = substring(names(MEs), 3) # name all the modules, from 3rd character on (first two are ME)
modNames## [1] "cyan" "black" "salmon" "pink" "red"
## [6] "green" "midnightblue" "blue" "brown" "purple"
## [11] "turquoise" "tan" "magenta" "lightcyan" "yellow"geneModuleMembership = as.data.frame(cor(vst2, MEs, use = "p"));
head(geneModuleMembership)## MEcyan MEblack MEsalmon
## Pocillopora_acuta_HIv2___RNAseq.g27841.t1 0.10438549 -0.12235265 -0.47423904
## Pocillopora_acuta_HIv2___RNAseq.g14011.t1 -0.16310820 0.21739302 -0.14153773
## Pocillopora_acuta_HIv2___RNAseq.g7479.t3a 0.08339074 -0.04725328 -0.33667736
## Pocillopora_acuta_HIv2___RNAseq.g7479.t3b 0.09142923 0.15329361 0.10385419
## Pocillopora_acuta_HIv2___RNAseq.g3340.t1 0.12842828 -0.10929165 0.07806305
## Pocillopora_acuta_HIv2___RNAseq.g682.t1 -0.02392754 -0.03296640 0.17552549
## MEpink MEred MEgreen
## Pocillopora_acuta_HIv2___RNAseq.g27841.t1 -0.3584060 -0.06200274 0.20465795
## Pocillopora_acuta_HIv2___RNAseq.g14011.t1 -0.2594287 -0.31902027 -0.01684585
## Pocillopora_acuta_HIv2___RNAseq.g7479.t3a 0.1186510 0.39184956 0.44854568
## Pocillopora_acuta_HIv2___RNAseq.g7479.t3b 0.5825581 0.82234053 0.35222463
## Pocillopora_acuta_HIv2___RNAseq.g3340.t1 0.2478598 0.35321199 0.60087070
## Pocillopora_acuta_HIv2___RNAseq.g682.t1 -0.4538716 -0.45358620 -0.31698652
## MEmidnightblue MEblue MEbrown
## Pocillopora_acuta_HIv2___RNAseq.g27841.t1 -0.21979419 0.49696801 0.2454424
## Pocillopora_acuta_HIv2___RNAseq.g14011.t1 0.08185613 0.08994886 -0.3855597
## Pocillopora_acuta_HIv2___RNAseq.g7479.t3a 0.14471353 0.48376147 0.3927873
## Pocillopora_acuta_HIv2___RNAseq.g7479.t3b 0.05633416 0.51529570 0.6432814
## Pocillopora_acuta_HIv2___RNAseq.g3340.t1 0.44897525 0.18646978 0.4855880
## Pocillopora_acuta_HIv2___RNAseq.g682.t1 -0.10512744 -0.50077743 -0.2318526
## MEpurple MEturquoise MEtan
## Pocillopora_acuta_HIv2___RNAseq.g27841.t1 0.35441847 -0.5622244 0.15275154
## Pocillopora_acuta_HIv2___RNAseq.g14011.t1 -0.41677235 0.1495377 -0.34664520
## Pocillopora_acuta_HIv2___RNAseq.g7479.t3a 0.34074711 -0.5934280 0.02119366
## Pocillopora_acuta_HIv2___RNAseq.g7479.t3b 0.33142341 -0.5855358 -0.20954032
## Pocillopora_acuta_HIv2___RNAseq.g3340.t1 0.39122310 -0.3553428 0.10450615
## Pocillopora_acuta_HIv2___RNAseq.g682.t1 -0.04828805 0.4289215 0.41719038
## MEmagenta MElightcyan MEyellow
## Pocillopora_acuta_HIv2___RNAseq.g27841.t1 0.07322422 0.28404105 0.53716214
## Pocillopora_acuta_HIv2___RNAseq.g14011.t1 0.16236103 0.15377483 0.01771476
## Pocillopora_acuta_HIv2___RNAseq.g7479.t3a -0.21754628 -0.06376141 0.33701433
## Pocillopora_acuta_HIv2___RNAseq.g7479.t3b -0.66545958 -0.41176150 0.08707934
## Pocillopora_acuta_HIv2___RNAseq.g3340.t1 -0.38686031 -0.34369924 -0.01298406
## Pocillopora_acuta_HIv2___RNAseq.g682.t1 0.26807670 0.31248662 0.01298851MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples));
names(geneModuleMembership) = paste("MM", modNames, sep="");
names(MMPvalue) = paste("p.MM", modNames, sep="");
# Flat treatment group
Flat = as.data.frame(Traits.Origin$Flat); # Define variable containing the desired column
names(Flat) = "Flat"
Flat_geneTraitSignificance = as.data.frame(cor(vst2, Flat, use = "p"));
Flat_GSPvalue = as.data.frame(corPvalueStudent(as.matrix(Flat_geneTraitSignificance), nSamples));
names(Flat_geneTraitSignificance) = paste("GS.", names(Flat), sep=""); # MA_geneTraitSignificance - pearsons correlation between reads and the MA grop
names(Flat_GSPvalue) = paste("p.GS.", names(Flat), sep=""); # corPvalueStudent
# Slope treatment group
Slope = as.data.frame(Traits.Origin$Slope); # Define variable containing the desired column
names(Slope) = "Slope"
Slope_geneTraitSignificance = as.data.frame(cor(vst2, Slope, use = "p"));
Slope_GSPvalue = as.data.frame(corPvalueStudent(as.matrix(Slope_geneTraitSignificance), nSamples));
names(Slope_geneTraitSignificance) = paste("GS.", names(Slope), sep=""); # MA_geneTraitSignificance - pearsons correlation between reads and the MA grop
names(Slope_GSPvalue) = paste("p.GS.", names(Slope), sep=""); # corPvalueStudentCall annotation data to get module gene data (prep for downstream GO)
library(rtracklayer)
gff <- rtracklayer::import("../../data/Pocillopora_acuta_HIv2.genes_fixed.gff3")
gff <- as.data.frame(gff)
dim(gff) # 478988 13## [1] 478988 13names(gff) ## [1] "seqnames" "start" "end" "width"
## [5] "strand" "source" "type" "score"
## [9] "phase" "ID" "transcript_id" "gene_id"
## [13] "Parent"transcripts <- subset(gff, type == "transcript")
transcripts_gr <- makeGRangesFromDataFrame(transcripts, keep.extra.columns=TRUE) #extract length information
transcript_lengths <- width(transcripts_gr) #isolate length of each gene
seqnames <- transcripts_gr$ID #extract list of gene id
lengths <- cbind(seqnames, transcript_lengths)
lengths <- as.data.frame(lengths) #convert to data frameeggnog <- read.delim("../../data/Pocillopora_acuta_HIv2.genes.EggNog_results.txt")
eggnog <- plyr::rename(eggnog, c("X.query"="gene_id"))gogene <- merge(transcripts, eggnog, by=c("gene_id"), all=T)colnames(gogene)## [1] "gene_id" "seqnames" "start" "end"
## [5] "width" "strand" "source" "type"
## [9] "score.x" "phase" "ID" "transcript_id"
## [13] "Parent" "seed_ortholog" "evalue" "score.y"
## [17] "eggNOG_OGs" "max_annot_lvl" "COG_category" "Description"
## [21] "Preferred_name" "GOs" "EC" "KEGG_ko"
## [25] "KEGG_Pathway" "KEGG_Module" "KEGG_Reaction" "KEGG_rclass"
## [29] "BRITE" "KEGG_TC" "CAZy" "BiGG_Reaction"
## [33] "PFAMs"probes = rownames(cts_filtpoa)
probes2annot = match(probes, gogene$gene_id)
# The following is the number or probes without annotation:
sum(is.na(probes2annot))## [1] 0# Should return 0.BUILD GENE INFO DATAFRAMES
length(probes)## [1] 9012length(gogene$seqnames[probes2annot])## [1] 9012length(gogene$seqnames[probes2annot])## [1] 9012# Create the starting data frame
names(Flat_geneTraitSignificance)## [1] "GS.Flat"names(Flat_GSPvalue)## [1] "p.GS.Flat"geneInfo_GROUPS = data.frame(substanceBXH = probes,
geneSymbol = gogene$seqnames[probes2annot],
moduleColor = moduleColors,
GO.terms = gogene$GOs[probes2annot],
GO.description = gogene$Description[probes2annot],
Flat_geneTraitSignificance, Slope_geneTraitSignificance, # call this specific to the module and trait of interest
Flat_GSPvalue, Slope_GSPvalue) # call this specific to the module and trait of interest
head(geneInfo_GROUPS)## substanceBXH
## Pocillopora_acuta_HIv2___RNAseq.g27841.t1 Pocillopora_acuta_HIv2___RNAseq.g27841.t1
## Pocillopora_acuta_HIv2___RNAseq.g14011.t1 Pocillopora_acuta_HIv2___RNAseq.g14011.t1
## Pocillopora_acuta_HIv2___RNAseq.g7479.t3a Pocillopora_acuta_HIv2___RNAseq.g7479.t3a
## Pocillopora_acuta_HIv2___RNAseq.g7479.t3b Pocillopora_acuta_HIv2___RNAseq.g7479.t3b
## Pocillopora_acuta_HIv2___RNAseq.g3340.t1 Pocillopora_acuta_HIv2___RNAseq.g3340.t1
## Pocillopora_acuta_HIv2___RNAseq.g682.t1 Pocillopora_acuta_HIv2___RNAseq.g682.t1
## geneSymbol
## Pocillopora_acuta_HIv2___RNAseq.g27841.t1 Pocillopora_acuta_HIv2___Sc0000000
## Pocillopora_acuta_HIv2___RNAseq.g14011.t1 Pocillopora_acuta_HIv2___Sc0000005
## Pocillopora_acuta_HIv2___RNAseq.g7479.t3a Pocillopora_acuta_HIv2___xfSc0000001
## Pocillopora_acuta_HIv2___RNAseq.g7479.t3b Pocillopora_acuta_HIv2___xfSc0000001
## Pocillopora_acuta_HIv2___RNAseq.g3340.t1 Pocillopora_acuta_HIv2___xfSc0000006
## Pocillopora_acuta_HIv2___RNAseq.g682.t1 Pocillopora_acuta_HIv2___xfSc0000002
## moduleColor
## Pocillopora_acuta_HIv2___RNAseq.g27841.t1 yellow
## Pocillopora_acuta_HIv2___RNAseq.g14011.t1 turquoise
## Pocillopora_acuta_HIv2___RNAseq.g7479.t3a blue
## Pocillopora_acuta_HIv2___RNAseq.g7479.t3b red
## Pocillopora_acuta_HIv2___RNAseq.g3340.t1 green
## Pocillopora_acuta_HIv2___RNAseq.g682.t1 turquoise
## GO.terms
## Pocillopora_acuta_HIv2___RNAseq.g27841.t1 GO:0000003,GO:0003006,GO:0003674,GO:0003824,GO:0004672,GO:0004674,GO:0004693,GO:0005575,GO:0005622,GO:0005623,GO:0005634,GO:0005737,GO:0005813,GO:0005815,GO:0005856,GO:0006464,GO:0006468,GO:0006793,GO:0006796,GO:0006807,GO:0007154,GO:0007165,GO:0007548,GO:0008150,GO:0008152,GO:0009987,GO:0015630,GO:0016301,GO:0016310,GO:0016740,GO:0016772,GO:0016773,GO:0019538,GO:0022414,GO:0023052,GO:0032502,GO:0036211,GO:0043170,GO:0043226,GO:0043227,GO:0043228,GO:0043229,GO:0043231,GO:0043232,GO:0043412,GO:0044237,GO:0044238,GO:0044260,GO:0044267,GO:0044422,GO:0044424,GO:0044430,GO:0044446,GO:0044464,GO:0050789,GO:0050794,GO:0050896,GO:0051716,GO:0051726,GO:0065007,GO:0071704,GO:0097472,GO:0140096,GO:1901564
## Pocillopora_acuta_HIv2___RNAseq.g14011.t1 GO:0000118,GO:0000123,GO:0003674,GO:0003676,GO:0003677,GO:0003712,GO:0003713,GO:0005488,GO:0005575,GO:0005622,GO:0005623,GO:0005634,GO:0005654,GO:0005730,GO:0006325,GO:0006338,GO:0006355,GO:0006357,GO:0006464,GO:0006473,GO:0006475,GO:0006807,GO:0006996,GO:0008150,GO:0008152,GO:0009889,GO:0009891,GO:0009893,GO:0009987,GO:0010468,GO:0010556,GO:0010557,GO:0010604,GO:0016043,GO:0016569,GO:0016570,GO:0016573,GO:0016604,GO:0016607,GO:0018193,GO:0018205,GO:0018393,GO:0018394,GO:0019219,GO:0019222,GO:0019538,GO:0030914,GO:0031248,GO:0031323,GO:0031325,GO:0031326,GO:0031328,GO:0031974,GO:0031981,GO:0032991,GO:0036211,GO:0043170,GO:0043226,GO:0043227,GO:0043228,GO:0043229,GO:0043231,GO:0043232,GO:0043233,GO:0043412,GO:0043543,GO:0043966,GO:0044237,GO:0044238,GO:0044260,GO:0044267,GO:0044422,GO:0044424,GO:0044428,GO:0044446,GO:0044451,GO:0044464,GO:0045935,GO:0048518,GO:0048522,GO:0050789,GO:0050794,GO:0051171,GO:0051173,GO:0051252,GO:0051254,GO:0051276,GO:0060255,GO:0065007,GO:0070013,GO:0070461,GO:0071704,GO:0071840,GO:0080090,GO:0097159,GO:0140110,GO:1901363,GO:1901564,GO:1902493,GO:1902494,GO:1902680,GO:1903506,GO:1903508,GO:1990234,GO:2000112,GO:2001141
## Pocillopora_acuta_HIv2___RNAseq.g7479.t3a <NA>
## Pocillopora_acuta_HIv2___RNAseq.g7479.t3b <NA>
## Pocillopora_acuta_HIv2___RNAseq.g3340.t1 GO:0003674,GO:0003824,GO:0005488,GO:0005515,GO:0005575,GO:0005622,GO:0005623,GO:0005737,GO:0005739,GO:0005740,GO:0005743,GO:0005759,GO:0006629,GO:0006720,GO:0006732,GO:0006733,GO:0006743,GO:0006744,GO:0006950,GO:0007275,GO:0007399,GO:0008150,GO:0008152,GO:0008299,GO:0008610,GO:0009058,GO:0009108,GO:0009987,GO:0010941,GO:0014016,GO:0014019,GO:0016020,GO:0016043,GO:0016740,GO:0016765,GO:0019866,GO:0019898,GO:0022008,GO:0022607,GO:0030154,GO:0031090,GO:0031312,GO:0031314,GO:0031966,GO:0031967,GO:0031974,GO:0031975,GO:0032501,GO:0032502,GO:0032991,GO:0033554,GO:0042180,GO:0042181,GO:0042981,GO:0043066,GO:0043067,GO:0043069,GO:0043226,GO:0043227,GO:0043229,GO:0043231,GO:0043233,GO:0043933,GO:0044085,GO:0044237,GO:0044238,GO:0044249,GO:0044255,GO:0044281,GO:0044283,GO:0044422,GO:0044424,GO:0044425,GO:0044429,GO:0044444,GO:0044446,GO:0044455,GO:0044464,GO:0046982,GO:0046983,GO:0048468,GO:0048519,GO:0048523,GO:0048699,GO:0048731,GO:0048856,GO:0048863,GO:0048864,GO:0048869,GO:0050789,GO:0050794,GO:0050896,GO:0051186,GO:0051188,GO:0051259,GO:0051262,GO:0051290,GO:0051291,GO:0051716,GO:0060548,GO:0065003,GO:0065007,GO:0070013,GO:0071704,GO:0071840,GO:0098780,GO:1901576,GO:1901661,GO:1901663,GO:1902494,GO:1990234
## Pocillopora_acuta_HIv2___RNAseq.g682.t1 <NA>
## GO.description
## Pocillopora_acuta_HIv2___RNAseq.g27841.t1 cyclin-dependent protein serine/threonine kinase activity
## Pocillopora_acuta_HIv2___RNAseq.g14011.t1 TAF6-like RNA polymerase II, p300 CBP-associated
## Pocillopora_acuta_HIv2___RNAseq.g7479.t3a <NA>
## Pocillopora_acuta_HIv2___RNAseq.g7479.t3b <NA>
## Pocillopora_acuta_HIv2___RNAseq.g3340.t1 trans-hexaprenyltranstransferase activity
## Pocillopora_acuta_HIv2___RNAseq.g682.t1 <NA>
## GS.Flat GS.Slope p.GS.Flat
## Pocillopora_acuta_HIv2___RNAseq.g27841.t1 0.2031928 -0.2031928 1.756178e-01
## Pocillopora_acuta_HIv2___RNAseq.g14011.t1 -0.3260066 0.3260066 2.703391e-02
## Pocillopora_acuta_HIv2___RNAseq.g7479.t3a 0.3491137 -0.3491137 1.740688e-02
## Pocillopora_acuta_HIv2___RNAseq.g7479.t3b 0.6161971 -0.6161971 5.140317e-06
## Pocillopora_acuta_HIv2___RNAseq.g3340.t1 0.3860958 -0.3860958 8.044440e-03
## Pocillopora_acuta_HIv2___RNAseq.g682.t1 -0.1301845 0.1301845 3.885043e-01
## p.GS.Slope
## Pocillopora_acuta_HIv2___RNAseq.g27841.t1 1.756178e-01
## Pocillopora_acuta_HIv2___RNAseq.g14011.t1 2.703391e-02
## Pocillopora_acuta_HIv2___RNAseq.g7479.t3a 1.740688e-02
## Pocillopora_acuta_HIv2___RNAseq.g7479.t3b 5.140317e-06
## Pocillopora_acuta_HIv2___RNAseq.g3340.t1 8.044440e-03
## Pocillopora_acuta_HIv2___RNAseq.g682.t1 3.885043e-01modOrder = order(-abs(cor(MEs, Flat, use = "p"))); # order by the strength of the correlation between module and trait values for each sample
for (mod in 1:ncol(geneModuleMembership)) # Add module membership information in the chosen order
{
oldNames = names(geneInfo_GROUPS)
geneInfo_GROUPS = data.frame(geneInfo_GROUPS, geneModuleMembership[, modOrder[mod]],
MMPvalue[, modOrder[mod]]);
names(geneInfo_GROUPS) = c(oldNames, paste("A.", modNames[modOrder[mod]], sep=""),
paste("p.A.", modNames[modOrder[mod]], sep=""))
}Order the genes in the geneInfo variable first by module color, then by geneTraitSignificance
geneOrder = order(geneInfo_GROUPS$moduleColor, -abs(geneInfo_GROUPS$GS.Flat));
geneInfo_GROUPS = geneInfo_GROUPS[geneOrder, ]
head(geneInfo_GROUPS)## substanceBXH
## Pocillopora_acuta_HIv2___RNAseq.g11459.t1 Pocillopora_acuta_HIv2___RNAseq.g11459.t1
## Pocillopora_acuta_HIv2___RNAseq.g27371.t1 Pocillopora_acuta_HIv2___RNAseq.g27371.t1
## Pocillopora_acuta_HIv2___RNAseq.30217_t Pocillopora_acuta_HIv2___RNAseq.30217_t
## Pocillopora_acuta_HIv2___RNAseq.g29237.t1 Pocillopora_acuta_HIv2___RNAseq.g29237.t1
## Pocillopora_acuta_HIv2___TS.g23133.t3 Pocillopora_acuta_HIv2___TS.g23133.t3
## Pocillopora_acuta_HIv2___RNAseq.g12199.t1 Pocillopora_acuta_HIv2___RNAseq.g12199.t1
## geneSymbol
## Pocillopora_acuta_HIv2___RNAseq.g11459.t1 Pocillopora_acuta_HIv2___Sc0000013
## Pocillopora_acuta_HIv2___RNAseq.g27371.t1 Pocillopora_acuta_HIv2___xfSc0000042
## Pocillopora_acuta_HIv2___RNAseq.30217_t Pocillopora_acuta_HIv2___xpSc0000388
## Pocillopora_acuta_HIv2___RNAseq.g29237.t1 Pocillopora_acuta_HIv2___xfSc0000022
## Pocillopora_acuta_HIv2___TS.g23133.t3 Pocillopora_acuta_HIv2___xfSc0000008
## Pocillopora_acuta_HIv2___RNAseq.g12199.t1 Pocillopora_acuta_HIv2___xfSc0000009
## moduleColor
## Pocillopora_acuta_HIv2___RNAseq.g11459.t1 black
## Pocillopora_acuta_HIv2___RNAseq.g27371.t1 black
## Pocillopora_acuta_HIv2___RNAseq.30217_t black
## Pocillopora_acuta_HIv2___RNAseq.g29237.t1 black
## Pocillopora_acuta_HIv2___TS.g23133.t3 black
## Pocillopora_acuta_HIv2___RNAseq.g12199.t1 black
## GO.terms
## Pocillopora_acuta_HIv2___RNAseq.g11459.t1 <NA>
## Pocillopora_acuta_HIv2___RNAseq.g27371.t1 <NA>
## Pocillopora_acuta_HIv2___RNAseq.30217_t -
## Pocillopora_acuta_HIv2___RNAseq.g29237.t1 <NA>
## Pocillopora_acuta_HIv2___TS.g23133.t3 GO:0000003,GO:0001101,GO:0001553,GO:0001868,GO:0001869,GO:0002020,GO:0002376,GO:0002437,GO:0002438,GO:0002526,GO:0002576,GO:0002673,GO:0002682,GO:0002683,GO:0002697,GO:0002698,GO:0002920,GO:0002921,GO:0003006,GO:0003674,GO:0004857,GO:0004866,GO:0004867,GO:0005096,GO:0005102,GO:0005488,GO:0005515,GO:0005575,GO:0005576,GO:0005615,GO:0005622,GO:0005623,GO:0005737,GO:0005829,GO:0006807,GO:0006810,GO:0006887,GO:0006950,GO:0006952,GO:0006953,GO:0006954,GO:0006955,GO:0007275,GO:0007548,GO:0007565,GO:0007566,GO:0007584,GO:0007596,GO:0007597,GO:0007599,GO:0008047,GO:0008150,GO:0008152,GO:0008406,GO:0008585,GO:0009605,GO:0009611,GO:0009719,GO:0009725,GO:0009790,GO:0009892,GO:0009966,GO:0009987,GO:0009991,GO:0010033,GO:0010037,GO:0010466,GO:0010468,GO:0010605,GO:0010629,GO:0010646,GO:0010951,GO:0010955,GO:0012505,GO:0014070,GO:0016043,GO:0016192,GO:0019222,GO:0019538,GO:0019838,GO:0019899,GO:0019955,GO:0019956,GO:0019958,GO:0019959,GO:0019966,GO:0022411,GO:0022414,GO:0022602,GO:0022607,GO:0022617,GO:0023051,GO:0030141,GO:0030154,GO:0030162,GO:0030198,GO:0030234,GO:0030414,GO:0030449,GO:0030695,GO:0031091,GO:0031093,GO:0031323,GO:0031324,GO:0031347,GO:0031348,GO:0031410,GO:0031667,GO:0031960,GO:0031974,GO:0031982,GO:0031983,GO:0032101,GO:0032268,GO:0032269,GO:0032501,GO:0032502,GO:0032504,GO:0032940,GO:0033993,GO:0034694,GO:0034695,GO:0034774,GO:0042060,GO:0042221,GO:0042493,GO:0042698,GO:0042802,GO:0042803,GO:0043062,GO:0043085,GO:0043086,GO:0043087,GO:0043120,GO:0043121,GO:0043170,GO:0043226,GO:0043227,GO:0043229,GO:0043233,GO:0043547,GO:0043933,GO:0044085,GO:0044092,GO:0044093,GO:0044238,GO:0044421,GO:0044422,GO:0044424,GO:0044433,GO:0044444,GO:0044446,GO:0044464,GO:0044703,GO:0044706,GO:0044877,GO:0045055,GO:0045088,GO:0045137,GO:0045824,GO:0045861,GO:0045916,GO:0046545,GO:0046660,GO:0046903,GO:0046983,GO:0048306,GO:0048403,GO:0048406,GO:0048511,GO:0048513,GO:0048519,GO:0048523,GO:0048545,GO:0048568,GO:0048583,GO:0048585,GO:0048608,GO:0048609,GO:0048731,GO:0048856,GO:0048863,GO:0048869,GO:0050727,GO:0050776,GO:0050777,GO:0050789,GO:0050790,GO:0050794,GO:0050817,GO:0050878,GO:0050896,GO:0051056,GO:0051171,GO:0051172,GO:0051179,GO:0051234,GO:0051246,GO:0051248,GO:0051259,GO:0051260,GO:0051262,GO:0051289,GO:0051336,GO:0051345,GO:0051346,GO:0051384,GO:0051704,GO:0052547,GO:0052548,GO:0060205,GO:0060255,GO:0060589,GO:0061134,GO:0061135,GO:0061458,GO:0065003,GO:0065007,GO:0065008,GO:0065009,GO:0070013,GO:0070542,GO:0070613,GO:0071704,GO:0071840,GO:0072347,GO:0072376,GO:0072378,GO:0080090,GO:0080134,GO:0097305,GO:0097708,GO:0098772,GO:0099503,GO:1901564,GO:1901654,GO:1901700,GO:1902531,GO:1903317,GO:1903318,GO:1990402,GO:2000257,GO:2000258
## Pocillopora_acuta_HIv2___RNAseq.g12199.t1 <NA>
## GO.description
## Pocillopora_acuta_HIv2___RNAseq.g11459.t1 <NA>
## Pocillopora_acuta_HIv2___RNAseq.g27371.t1 <NA>
## Pocillopora_acuta_HIv2___RNAseq.30217_t GTPase activator activity
## Pocillopora_acuta_HIv2___RNAseq.g29237.t1 <NA>
## Pocillopora_acuta_HIv2___TS.g23133.t3 endopeptidase inhibitor activity
## Pocillopora_acuta_HIv2___RNAseq.g12199.t1 <NA>
## GS.Flat GS.Slope p.GS.Flat
## Pocillopora_acuta_HIv2___RNAseq.g11459.t1 0.5677405 -0.5677405 3.872268e-05
## Pocillopora_acuta_HIv2___RNAseq.g27371.t1 0.5452501 -0.5452501 8.913481e-05
## Pocillopora_acuta_HIv2___RNAseq.30217_t 0.5127223 -0.5127223 2.693006e-04
## Pocillopora_acuta_HIv2___RNAseq.g29237.t1 0.4741063 -0.4741063 8.732422e-04
## Pocillopora_acuta_HIv2___TS.g23133.t3 0.4645548 -0.4645548 1.144263e-03
## Pocillopora_acuta_HIv2___RNAseq.g12199.t1 0.4551104 -0.4551104 1.483628e-03
## p.GS.Slope A.brown p.A.brown
## Pocillopora_acuta_HIv2___RNAseq.g11459.t1 3.872268e-05 0.3735665 0.01055044
## Pocillopora_acuta_HIv2___RNAseq.g27371.t1 8.913481e-05 0.3655440 0.01248495
## Pocillopora_acuta_HIv2___RNAseq.30217_t 2.693006e-04 0.2255392 0.13179089
## Pocillopora_acuta_HIv2___RNAseq.g29237.t1 8.732422e-04 0.3931221 0.00687840
## Pocillopora_acuta_HIv2___TS.g23133.t3 1.144263e-03 0.3318744 0.02424532
## Pocillopora_acuta_HIv2___RNAseq.g12199.t1 1.483628e-03 0.2863948 0.05365629
## A.magenta p.A.magenta A.red
## Pocillopora_acuta_HIv2___RNAseq.g11459.t1 -0.4390326 2.270613e-03 0.3067458
## Pocillopora_acuta_HIv2___RNAseq.g27371.t1 -0.5153303 2.474598e-04 0.2811921
## Pocillopora_acuta_HIv2___RNAseq.30217_t -0.4587281 1.344296e-03 0.4699351
## Pocillopora_acuta_HIv2___RNAseq.g29237.t1 -0.6079115 7.434349e-06 0.5714810
## Pocillopora_acuta_HIv2___TS.g23133.t3 -0.5510449 7.231490e-05 0.3696800
## Pocillopora_acuta_HIv2___RNAseq.g12199.t1 -0.5088948 3.045180e-04 0.4181299
## p.A.red A.turquoise
## Pocillopora_acuta_HIv2___RNAseq.g11459.t1 3.812803e-02 -0.10583159
## Pocillopora_acuta_HIv2___RNAseq.g27371.t1 5.835180e-02 -0.21881675
## Pocillopora_acuta_HIv2___RNAseq.30217_t 9.835954e-04 -0.09531465
## Pocillopora_acuta_HIv2___RNAseq.g29237.t1 3.350908e-05 -0.04948723
## Pocillopora_acuta_HIv2___TS.g23133.t3 1.145282e-02 -0.23930488
## Pocillopora_acuta_HIv2___RNAseq.g12199.t1 3.832660e-03 -0.03439068
## p.A.turquoise A.purple p.A.purple
## Pocillopora_acuta_HIv2___RNAseq.g11459.t1 0.4839269 -0.001344154 0.9929264
## Pocillopora_acuta_HIv2___RNAseq.g27371.t1 0.1440109 0.097484918 0.5192403
## Pocillopora_acuta_HIv2___RNAseq.30217_t 0.5286280 -0.242928037 0.1037909
## Pocillopora_acuta_HIv2___RNAseq.g29237.t1 0.7439690 -0.085020600 0.5742610
## Pocillopora_acuta_HIv2___TS.g23133.t3 0.1092114 -0.051265727 0.7350959
## Pocillopora_acuta_HIv2___RNAseq.g12199.t1 0.8205039 -0.060136362 0.6913677
## A.green p.A.green A.lightcyan
## Pocillopora_acuta_HIv2___RNAseq.g11459.t1 0.23766284 1.117381e-01 -0.1952640
## Pocillopora_acuta_HIv2___RNAseq.g27371.t1 0.26820482 7.151857e-02 -0.1217634
## Pocillopora_acuta_HIv2___RNAseq.30217_t 0.16609603 2.699402e-01 -0.2728921
## Pocillopora_acuta_HIv2___RNAseq.g29237.t1 -0.04119175 7.857737e-01 -0.3495440
## Pocillopora_acuta_HIv2___TS.g23133.t3 0.55430775 6.417267e-05 -0.3123909
## Pocillopora_acuta_HIv2___RNAseq.g12199.t1 0.20017449 1.822589e-01 -0.4310812
## p.A.lightcyan A.pink p.A.pink
## Pocillopora_acuta_HIv2___RNAseq.g11459.t1 0.193444104 0.18853025 0.2095613030
## Pocillopora_acuta_HIv2___RNAseq.g27371.t1 0.420174721 0.08986363 0.5525688928
## Pocillopora_acuta_HIv2___RNAseq.30217_t 0.066519359 0.37522665 0.0101840079
## Pocillopora_acuta_HIv2___RNAseq.g29237.t1 0.017259628 0.53290057 0.0001374455
## Pocillopora_acuta_HIv2___TS.g23133.t3 0.034545563 0.25217327 0.0908943212
## Pocillopora_acuta_HIv2___RNAseq.g12199.t1 0.002781693 0.44616937 0.0018845316
## A.blue p.A.blue A.salmon
## Pocillopora_acuta_HIv2___RNAseq.g11459.t1 0.03781604 0.8029651 0.4435399
## Pocillopora_acuta_HIv2___RNAseq.g27371.t1 0.16078732 0.2857617 0.3429064
## Pocillopora_acuta_HIv2___RNAseq.30217_t 0.18424167 0.2203025 0.3152706
## Pocillopora_acuta_HIv2___RNAseq.g29237.t1 0.07874465 0.6029359 0.4624926
## Pocillopora_acuta_HIv2___TS.g23133.t3 0.27193253 0.0675195 0.2510276
## Pocillopora_acuta_HIv2___RNAseq.g12199.t1 -0.03636981 0.8103591 0.4161561
## p.A.salmon A.midnightblue
## Pocillopora_acuta_HIv2___RNAseq.g11459.t1 0.002019420 0.34284439
## Pocillopora_acuta_HIv2___RNAseq.g27371.t1 0.019652076 0.22524824
## Pocillopora_acuta_HIv2___RNAseq.30217_t 0.032827874 0.24143271
## Pocillopora_acuta_HIv2___RNAseq.g29237.t1 0.001211796 0.01104787
## Pocillopora_acuta_HIv2___TS.g23133.t3 0.092421359 0.55037413
## Pocillopora_acuta_HIv2___RNAseq.g12199.t1 0.004020175 0.33515410
## p.A.midnightblue A.black
## Pocillopora_acuta_HIv2___RNAseq.g11459.t1 1.967570e-02 0.5029437
## Pocillopora_acuta_HIv2___RNAseq.g27371.t1 1.323031e-01 0.4355470
## Pocillopora_acuta_HIv2___RNAseq.30217_t 1.060024e-01 0.6425922
## Pocillopora_acuta_HIv2___RNAseq.g29237.t1 9.419093e-01 0.4207111
## Pocillopora_acuta_HIv2___TS.g23133.t3 7.410146e-05 0.4375792
## Pocillopora_acuta_HIv2___RNAseq.g12199.t1 2.279446e-02 0.4473112
## p.A.black A.cyan p.A.cyan
## Pocillopora_acuta_HIv2___RNAseq.g11459.t1 3.675918e-04 0.200348223 0.18187183
## Pocillopora_acuta_HIv2___RNAseq.g27371.t1 2.483421e-03 0.273605353 0.06578363
## Pocillopora_acuta_HIv2___RNAseq.30217_t 1.474083e-06 -0.009882091 0.94803020
## Pocillopora_acuta_HIv2___RNAseq.g29237.t1 3.599055e-03 -0.045947451 0.76172625
## Pocillopora_acuta_HIv2___TS.g23133.t3 2.357306e-03 0.179255763 0.23326230
## Pocillopora_acuta_HIv2___RNAseq.g12199.t1 1.828487e-03 0.149449675 0.32153753
## A.yellow p.A.yellow A.tan
## Pocillopora_acuta_HIv2___RNAseq.g11459.t1 -0.04482242 0.76739590 -0.0465158
## Pocillopora_acuta_HIv2___RNAseq.g27371.t1 0.13687184 0.36438494 -0.1011951
## Pocillopora_acuta_HIv2___RNAseq.30217_t -0.10080009 0.50506242 -0.4158460
## Pocillopora_acuta_HIv2___RNAseq.g29237.t1 -0.21492986 0.15144861 -0.3403418
## Pocillopora_acuta_HIv2___TS.g23133.t3 0.03639833 0.81021311 -0.3617676
## Pocillopora_acuta_HIv2___RNAseq.g12199.t1 -0.27200931 0.06743904 -0.2045885
## p.A.tan
## Pocillopora_acuta_HIv2___RNAseq.g11459.t1 0.758866700
## Pocillopora_acuta_HIv2___RNAseq.g27371.t1 0.503386444
## Pocillopora_acuta_HIv2___RNAseq.30217_t 0.004050348
## Pocillopora_acuta_HIv2___RNAseq.g29237.t1 0.020648495
## Pocillopora_acuta_HIv2___TS.g23133.t3 0.013495757
## Pocillopora_acuta_HIv2___RNAseq.g12199.t1 0.172606503write.csv(geneInfo_GROUPS, file = "../../output/WGCNA/WGCNA_ModuleMembership.csv")Individual Module Plots again, with n = number of genes in module
nrow(geneInfo_GROUPS %>% dplyr::filter(moduleColor=="brown")) #942## [1] 942nrow(geneInfo_GROUPS %>% dplyr::filter(moduleColor=="red")) #425## [1] 425nrow(geneInfo_GROUPS %>% dplyr::filter(moduleColor=="turquoise")) #2558## [1] 2558nrow(geneInfo_GROUPS %>% dplyr::filter(moduleColor=="magenta")) #219## [1] 219Brown only - By Treatment and Origin, Plot mean module eigengene
n_mod <- nrow(geneInfo_GROUPS %>% dplyr::filter(moduleColor=="brown")) #942
pdf("../../output/WGCNA/plot_MEs_fig_Trt_Orig_BROWN.pdf")
plot_MEs_fig_Trt_Orig <- ggplot(filter(plot_MEs, Module =="MEbrown"), aes(y=Mean, x=Treatment, color=Origin, fill=Origin))+
geom_boxplot(alpha=0.4, color="black", outlier.shape = NA, size =1)+
geom_point(alpha=0.8,position=position_jitterdodge(0.05),size=2.5)+
geom_hline(yintercept = 0, linetype="dashed", color = 'black', size=0.7, show.legend = TRUE)+
scale_fill_manual("Origin", values=c("Slope"='paleturquoise3', "Flat"= "indianred"))+
scale_color_manual("Origin", values=c("Slope"='paleturquoise3', "Flat"= "indianred"))+
scale_y_continuous(expression(Mean~module~eigenegene), limits=c(-0.4, 0.4), expand = c(0.005, 0.005))+
theme_classic()+
labs(title = "Brown Module") +
annotate("text", x = Inf, y = Inf, label = paste("n genes =", n_mod), hjust = 1.1, vjust = 2.1, size = 6, color = "black") +
theme(axis.text.x=element_text(vjust=0.5, size=20),
axis.text.y=element_text(vjust=0.5, size=20),
axis.title.x=element_blank(),
axis.title.y=element_text(size=20),
legend.text = element_text(vjust=0.5, size=20),
panel.background= element_rect(fill=NA, color='black'),
strip.text = element_text(vjust=0.5, size=20),
legend.position = "none")
plot_MEs_fig_Trt_Orig
dev.off()## quartz_off_screen
## 2plot_MEs_fig_Trt_Orig
Red only - By Treatment and Origin, Plot mean module eigengene
n_mod <- nrow(geneInfo_GROUPS %>% dplyr::filter(moduleColor=="red")) #425
pdf("../../output/WGCNA/plot_MEs_fig_Trt_Orig_RED.pdf")
plot_MEs_fig_Trt_Orig <- ggplot(filter(plot_MEs, Module =="MEred"), aes(y=Mean, x=Treatment, color=Origin, fill=Origin))+
geom_boxplot(alpha=0.4, color="black", outlier.shape = NA, size =1)+
geom_point(alpha=0.8,position=position_jitterdodge(0.05),size=2.5)+
geom_hline(yintercept = 0, linetype="dashed", color = 'black', size=0.7, show.legend = TRUE)+
scale_fill_manual("Origin", values=c("Slope"='paleturquoise3', "Flat"= "indianred"))+
scale_color_manual("Origin", values=c("Slope"='paleturquoise3', "Flat"= "indianred"))+
scale_y_continuous(expression(Mean~module~eigenegene), limits=c(-0.4, 0.4), expand = c(0.005, 0.005))+
theme_classic()+
labs(title = "Red Module")+
annotate("text", x = Inf, y = Inf, label = paste("n genes =", n_mod), hjust = 1.1, vjust = 2.1, size = 6, color = "black") +
theme(axis.text.x=element_text(vjust=0.5, size=20),
axis.text.y=element_text(vjust=0.5, size=20),
axis.title.x=element_blank(),
axis.title.y=element_text(size=20),
legend.text = element_text(vjust=0.5, size=20),
panel.background= element_rect(fill=NA, color='black'),
strip.text = element_text(vjust=0.5, size=20),
legend.position = "none")
plot_MEs_fig_Trt_Orig
dev.off()## quartz_off_screen
## 2plot_MEs_fig_Trt_Orig
Turquoise only - By Treatment and Origin, Plot mean module eigengene
n_mod <- nrow(geneInfo_GROUPS %>% dplyr::filter(moduleColor=="turquoise")) #2558
pdf("../../output/WGCNA/plot_MEs_fig_Trt_Orig_TURQUOISE.pdf")
plot_MEs_fig_Trt_Orig <- ggplot(filter(plot_MEs, Module =="MEturquoise"), aes(y=Mean, x=Treatment, color=Origin, fill=Origin))+
geom_boxplot(alpha=0.4, color="black", outlier.shape = NA, size =1)+
geom_point(alpha=0.8,position=position_jitterdodge(0.05),size=2.5)+
geom_hline(yintercept = 0, linetype="dashed", color = 'black', size=0.7, show.legend = TRUE)+
scale_fill_manual("Origin", values=c("Slope"='paleturquoise3', "Flat"= "indianred"))+
scale_color_manual("Origin", values=c("Slope"='paleturquoise3', "Flat"= "indianred"))+
scale_y_continuous(expression(Mean~module~eigenegene), limits=c(-0.4, 0.4), expand = c(0.005, 0.005))+
theme_classic()+
labs(title = "Turquoise Module")+
annotate("text", x = Inf, y = Inf, label = paste("n genes =", n_mod), hjust = 1.1, vjust = 2.1, size = 6, color = "black") +
theme(axis.text.x=element_text(vjust=0.5, size=20),
axis.text.y=element_text(vjust=0.5, size=20),
axis.title.x=element_blank(),
axis.title.y=element_text(size=20),
legend.text = element_text(vjust=0.5, size=20),
panel.background= element_rect(fill=NA, color='black'),
strip.text = element_text(vjust=0.5, size=20),
legend.position = "none")
plot_MEs_fig_Trt_Orig
dev.off()## quartz_off_screen
## 2plot_MEs_fig_Trt_Orig
Magenta only - By Treatment and Origin, Plot mean module eigengene
n_mod <- nrow(geneInfo_GROUPS %>% dplyr::filter(moduleColor=="magenta")) #219
pdf("../../output/WGCNA/plot_MEs_fig_Trt_Orig_MAGENTA.pdf")
plot_MEs_fig_Trt_Orig <- ggplot(filter(plot_MEs, Module =="MEmagenta"), aes(y=Mean, x=Treatment, color=Origin, fill=Origin))+
geom_boxplot(alpha=0.4, color="black", outlier.shape = NA, size =1)+
geom_point(alpha=0.8,position=position_jitterdodge(0.05),size=2.5)+
geom_hline(yintercept = 0, linetype="dashed", color = 'black', size=0.7, show.legend = TRUE)+
scale_fill_manual("Origin", values=c("Slope"='paleturquoise3', "Flat"= "indianred"))+
scale_color_manual("Origin", values=c("Slope"='paleturquoise3', "Flat"= "indianred"))+
scale_y_continuous(expression(Mean~module~eigenegene), limits=c(-0.4, 0.4), expand = c(0.005, 0.005))+
theme_classic()+
labs(title = "Magenta Module")+
annotate("text", x = Inf, y = Inf, label = paste("n genes =", n_mod), hjust = 1.1, vjust = 2.1, size = 6, color = "black") +
theme(axis.text.x=element_text(vjust=0.5, size=20),
axis.text.y=element_text(vjust=0.5, size=20),
axis.title.x=element_blank(),
axis.title.y=element_text(size=20),
legend.text = element_text(vjust=0.5, size=20),
panel.background= element_rect(fill=NA, color='black'),
strip.text = element_text(vjust=0.5, size=20),
legend.position = "none")
plot_MEs_fig_Trt_Orig
dev.off()## quartz_off_screen
## 2plot_MEs_fig_Trt_Orig